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I.  INIFODOCTION 


Since  World  War  II,  combat  modeling,  simulation  and 
analysis  have  been  the  subjects  of  considerable  research. 
The  objectives  of  this  research  are  to  support  defense  deci¬ 
sion  making  and  doctrinal  developments  during  peace  and  war 
time.  During  peacetime  defense-planners  are  primarily 
concerned  with  weapon  procurement,  development,  acquisition, 
organisation  and  structuring.  During  war  time  it  is 
believed  that  a  better  understanding  of  the  quantitative 
aspects  cf  attrition  can  help  commanders  make  better  command 
and  control  decisions. 

Combat  processes  involve  complicated  interactions 
between  opposing  forces.  These  interactions  are  often  influ¬ 
enced  by  many  external  factors  such  as  environment,  troop 
quality  and  tactics.  There  are  different  types  of  combat 
models  such  as  war  games,  simulations  and  analytical  models. 
A  fundamental  requirement  for  a  good  model  is  that  it  must 
be  of  a  fairly  high  degree  of  operational  realism,  since 
otherwise  they  would  not  be  credible  to  military  planners. 
Cn  the  other  hand,  excessively  complicated  models  can  make 
the  mathematics  too  difficult  tc  handle. 

In  this  thesis,  a  generalised  Lanchester  [Hef.  1]  model 
which  contains  area-fire,  aimed-fire,  self-attrition  and 
replenishment  coefficients  is  used.  It  consists  of  a  system 
of  2N  bilinear  eguations  and  belongs  to  the  general  category 
of  analytical  models.  The  model  is  rich  enough  to  treat 
modern  combined-arms  operations  involving  heterogeneous 
forces.  It  is  also  possible  to  extend  the  model  to  analyse 
operations  cn  two  or  more  fronts. 

Among  the  many  important  issues  that  could  be  analysed 
using  this  model,  the  problem  of  optimum  force  distribution 


had  teen  studied  by  Bozencraf  t  and  Moose  (1983)  .  In  their 
paper  [Ref.  2],  an  objective  function  was  chosen  as  the 
difference  of  the  aggregate  attrition  rates.  It  was  shewn 
that  the  optimization  problem  is  mathematically  equivalent 
to  a  matrix  game.  Hence,  the  model  has  a  saddle-point  solu¬ 
tion  with  corresponding  optimum  force  distribution  vectors 
x  and  y  for  Blue  and  Crange  forces  respectively. 

In  addition,  the  neighborhood  stability  of  the  model  at 
the  operating  point  (x*  and  y*)  was  also  investigated.  By 
defining  two  parameters,  K1  and  K2  which  are  obtained  by 
considering  small  perturbations  around  the  operating  points, 
a  great  deal  could  be  learned  about  stability. 

activated  by  these  results,  much  of  the  work  done  during 
the  initial  part  of  this  thesis  was  directed  at  studying  the 
effect  of  stability  cn  battle  outcome.  The  ultimate  guestion 
is,  how  do  we  exploit  the  knowledge  of  stability  of  an  oper¬ 
ating  point  to  influence  battle  outcome?  Before  this  gues¬ 
tion  can  he  answered,  it  appears  that  there  is  a  need  for  a 
tetter  understanding  of  the  equilibrium  points.  Chapter  III 
is  devoted  to  finding  and  understanding  the  equilibrium 
solutions  and  their  stability  behavior.  Like  many  other 
nonlinear  system  of  equations,  the  Lanchester's  model 
adopted  here  has  multiple  equilibria.  Stability  analysis 
[Bef.  3]  of  a  non-linear  system  is  usually  done  by  methods 
which  do  net  require  prior  knowledge  of  the  equilibrium 
solutions.  One  example  of  such  a  method  is  the  Liapunov 
method  [Ref.  4],  If,  by  some  realizable  means,  the  equilib¬ 
rium  solutions  can  be  found  explicitly  then  there  is  no  need 
to  rely  on  these  indirect  methods  which  are  often  difficult 
to  implement. 

One  of  the  reasons  for  resorting  to  the  Liapunov  method 
is  the  difficulty  in  obtaining  equilibrium  solutions  of  a 
non-linear  system.  Many  numerical  methods  are  unsuitable  for 
reasons  such  as  difficulty  in  obtaining  good  initial 


13 


guesses,  non-convergence,  ill-conditioning  and  so  forth. 
Fortunately,  a  powerful  numerical  technigue  called  the 
Continuation  Method  can  be  applied  for  our  purpose.  This 
method  not  only  finds  all  the  solutions  (i.e.  it  is  exhaus¬ 
tive),  it  dees  not  even  require  initial  guesses. 

In  order  to  gain  a  firm  grasp  on  the  dynamics  of  the 
system  surrounding  the  equilibria,  it  is  helpful  tc  tempo¬ 
rarily  focus  attention  on  the  homogeneous  (1*1)  system.  In 
spite  of  its  simplicity,  the  1*1  system  is  not  devoid  of  the 
essential  characteristics  of  the  N*N  system.  In  fact,  the 
1*1  model  is  sufficiently  sophisticated  for  certain  analyses 
in  which  the  opposing  forces  can  be  assumed  to  be  homoge¬ 
neous.  As  we  proceed  through  Chapter  IV,  it  will  become 
clear  that  much  insight  into  the  stability  and  system 
dynamics  could  be  gained  by  merely  considering  the  1*1 
system.  Part  of  the  chapter  is  devoted  to  the  derivations 
and  interpretations  of  the  relations  between  system  asymp¬ 
totes,  locations  of  equilibrium  points  and  stability.  The 
dynamics  of  the  system  are  studied  using  the  idea  of  phase 
trajectories.  These  trajectories  represent  changes  of  force 
levels  with  time  and  they  will  be  shown  to  depend  not  only 
on  the  stabilities  of  equilibrium  points  but  also  on  the 
domains  of  attraction. 

Chapter  V  concentrates  on  tattle  outcome  which  is  one  of 
the  main  issues  facing  a  commander.  It  encompasses  many 
issues  such  as,  (1)  Who  will  win  and  by  what  margin?  (2) 
What  is  the  length  of  battle?  (3)  How  do  initial  deployments 
affect  battle  outcome?  (4)  Which  parameters  affect  battle 
outcome  most?  But  we  will  only  address  the  two  following 
subjects  : 

(a)  The  effect  of  stability  on  battle  outcome; 

(b)  The  effect  cf  varying  X  and  Y,  the  initial  force 

levels . 


The  lasic  approach  is  to  define  a  multistage  battle  with 
a  predetermined  condition  for  termination.  The  resultant 
payoff  matrix  can  then  be  used  to  obtain  the  optimum  set  of 
mixed  strategies.  An  example,  which  employs  KOREAN  WAR  data, 
is  presented  for  the  purpose  of  illustrations  and 
discussions. 

The  essence  of  the  findings  are: 

1.  Unstable  operating  conditions  can  be  exploited  to 
influence  battle  outcome,  especially  when  total  war 
resources  are  large.  The  effect  on  battle  outcome  is 
mere  pronounced  for  highly  unstable  warfare; 

2.  Initial  force  deployment  can  be  optimized  in  accor¬ 
dance  with  a  set  of  mixed  strategies. 

We  conclude  this  introduction  by  stating  two  of  the 
outstanding  issues.  The  first  question  is  the  extent  to 
which  one  can  replace  the  N*N  problem  by  the  1*1  problem. 
The  motivation  to  find  an  equivalent  1*1  system  stems  from 
(1)  our  better  understanding  of  the  1*1  system,  (2)  ease  of 
presentieg  and  visualizing  two-dimensional  pictures,  and  (3) 
savings  in  computational  effort. 

The  second  question  concerns  replenishment  rates.  In 
this  thesis,  the  replenishment  terms  used  in  the  model  have 
been  constant.  It  is  therefore  reasonable  to  ask,  how  to 
modify  replenishment  terms  to  reflect  a  higher  degree  of 
operational  realism?  In  other  words,  are  there  more  suit¬ 
able  time-dependent  replenishment  rates  r(t)? 
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II.  IANCHE STEE1 S  EQUATION 


3 

i 


A.  EACKGBCUND 

Combat  models  have  been  studied  as  a  form  of  decision 
aid  for  defense  planting.  A  wide  variety  of  defense  plan¬ 
ning  problems,  ranging  from  force  structuring  and  weapon 
selection  to  rates  of  deployment  in  battles  have  been  anal¬ 
ysed  vising  combat  models.  There  are  many  different  types  of 
models.  They  can  be  loosely  categorized  as  either  war  games, 
simulations  or  analytical  models.  Discussions  on  the 
nature,  advantages  and  shortcomings  of  each  can  be  found  in 
[Bef.  5  j. 

Cur  attention  will  be  focused  on  a  generalized 

lanchester's  [Ref.  5]  model,  which  is  an  analytical  model. 
It  consists  basically  of  a  system  of  ordinary  differential 
equations  describing  the  mutual  interactions  between 
opposing  combat  forces.  Although  earlier  works  in 

lanchester’s  model  [Bef.  6]  employed  only  a  few  terms  in  the 
equations,  modern  high  speed  computers  enable  more  general¬ 
ised,  realistic  and  responsive  versions  to  be  used. 

Consider  a  battlefield  with  opposing  forces,  Elue  and 
Crange,  denoted  by  { }  and  {  y ^ }  respectively.  The 

subscripts  i,  j  refer  to  the  type  of  forces  such  as 

infantry,  tanks,  artillery,  etc.  A  generalised  version  of 
lanchester's  model  given  by 


=  -xlUi  -  -  £bij*j  *  ri 


y.  =  -v.y.  -  y.y^x.c.  .  -  •  d .  .  +  s. 

J  j  '  li]  Li  l  i]  3 


( eqn  2.1) 


1  ^  I 

I  *>  T 

1  >-»•••* 


where 


ui  ’  vj 

a  .  .  ,  c  .  . 
13  13 

bir  dj. 

ri  ■  sj 


=  self-attrition  coefficients 
=  area-fire  attrition  coefficients 
r  aimed-fire  attrition  coefficients 
=  replenishment  coefficients 


is  adopted  in  this  thesis. 

Note  that  in  general  I  ?  Jf  implying  that  the  force 
compositions  may  be  different  for  the  two  sides.  It  is  also 
possible  to  extend  the  above  formulation  to  a  scenerio 
involving  more  than  cne  battlefield. 

In  the  next  two  sections,  the  highlights  of  the  work 
done  by  Hozencraft  and  Moose  (1983)  are  given.  The  work  done 
in  this  thesis  is  a  continuation  and  extention  of  their 
work.  The  detailed  derivations  of  the  results  obtained  by 
them  can  be  found  in  [ Hef .  2],  and  hence  are  not  included 
here. 

B.  OPTIMUM  FORCE  DISTRIBUTION 

The  question  of  optimum  force  distribution  arises  in 
combined-arms  operations.  The  problem  is  fundamentally 
this:  Given  aggregated  forces  X,  y,  how  should  one 
distribute  them  among  the  different  types  x^^  and  y^  ,  .i  = 
1,2..., I,  3  =  1,2...,J?  Since  loss  rate  is  one  of  the 
fundamental  concepts  in  combat  modeling,  it  is  reasonable  to 
choose  this  measure  as  a  starting  point.  The  objective 
function  was  chosen  to  be 
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M  i  -  ri’  -  Jj?i  -  sj 


(ecn  2.2) 


Fcr  this  choice  of  M,  it  was  shown  that  there  exists 
optimum  force  distribution  (row  and  column)  vectors  x*  and 
y*  such  that  for  any  ether  vectors  x  and  y 


~  *  *  *  A 

xAy  <  M  <  x  Ay 


(eqn  2.3) 


where 


*  *  * 

x  Av 


A  =  matrix  determined  by  attrition  coefficients  and 
the  aggregate  force  levels  X  and  V 

Ihe  resemblance  of  this  result  to  the  Mini  max  theorem 
[Hef.  7]  in  matrix  games  is  very  striking.  Indeed,  this 
result  holds  precisely  because  w  can  be  written  in  a  form 
mathematically  eguivalent  tc  a  matrix  game.  consequently, 
it  is  net  surprising  that  one  can  solve  for  the  optimum 
vector  x*  and  y*  by  means  of  a  linear  Program.  An  interac¬ 
tive  program  to  solve  a  2*2  program  is  given  in  Appendix  A. 

C.  BEIGEBCBHOOD  STAEI1ITY 

Equilibrium  conditions  can  be  achieved  if  the  replenish¬ 
ment  rates  are  chosen  to  make 


at  x  =  x*  and  y  =  y* .  Following  the  usual  approach  in  the 
analysis  of  nonlinear  system  stability,  equation  2.1  can 
then  be  transformed  into  a  system  of  linear  equations. 


<5x 

<5x 

i — 

<  CQ 

<  < 

1 _ 

• 

=  -C 

>  C  = 

A  A 

6y 

<5v 

D  Cl 

C  is  called  the  conflict  matrix  and  its  elements  are  deter¬ 
mined  by  the  attrition  coefficients  and  the  optimum  vectors 

A  A 

x*  and  y* .  For  the  system  of  equations  2.1,  A  and  C  are 
diagonal  matrices.  It  was  shown  that  two  parameters  kx,  k2 
partially  characterize  the  stability  of  the  system.  k1  and  k 
turn  out  to  be  the  column  sums  of  the  left  and  right  side  o 
the  matrix 


Eenoting  the  elements  of  the  sutmatrices  A,  B,  c,  D,  by  a^, 
bXj  ,  Cjj  and  dli  respectively,  kL  and  k2  can  be  written  as 


1 


independent  of  the  columns  i,  j.  Furthermore,  it  was  shown 
that  the  following  relation  holds 

5X  -  k.6X  =  6  Y  -  k-,6Y 


19 


Ml  K) 


It  was  found  that  the  equilibrium  point  (x*,  y*  }  is 
stable  if  kx  and  k2  are  negative.  If  k1  and  k2  are  positive, 
then  the  system  is  ’unstable'.1  Furthermore,  values  cf 
and  k2  and  hence  the  stability  of  the  operating  point  was 
found  to  be  affected  ty  the  aggregate  X  and  X. 


»f!ore  generally,  it  can  be.shown  that  k-,  <  ,  where  ' 

is  the  maximum  eigenvalue  of  -C. 


III.  MOIUDIME  NSION AL  (N * N)  SYSTEM 
i.  H1TOBE  CF  N*  H  PRCELEM 

The  interesting  results  highlighted  in  the  last  chapter 
provided  motivation  tc  extend  the  body  of  knowledge.  A  study 
of  the  effect  on  stability  of  battle  outcome  seems  to  have 
important  potentials  for  applications.  Should  a  commander 
strive  to  establish  a  stable  operating  point,  and  if  so, 
under  what  conditions?  Also,  what  is  the  optimum  initial 
level  of  forces  he  should  deploy  and  how  many  should  he 
maintain  in  reserve?  To  answer  these  questions,  more  knowl¬ 
edge  about  the  nature  of  these  equilibria  and  their 
stability  behavior  is  required. 

The  next  section  outlines  the  kind  of  problems  we  would 
expect  to  see  and  their  potential  complexity.  It  is 
followed  by  a  section  on  finding  the  equilibrium  solutions. 

1  •  Existence  of  Eultip le  Equilibria 

An  N*N  system  is  in  equilibrium  if  the  replenishment 

rates  ,  s.  are  such  that  there  is  no  change  in  the  force 

•  ^  • 

levels  (xi  =  y^  =  0) .  The  system  of  equations  becomes 


0  =  -x-u-  -  x-Va-.y.  -  Vb-  +  r. 

ii  i  i  j  '  j  ^  j  l 

j  j 

0  *  -  v  •  y  .  -  y  •  T'x  •  c  .  -  *  Vx  .  d  .  .  +  s  . 

JT  J  }  J  ij  A.  i  ij  j 


( eqn  3.1) 


i,j  =  1,2, 


N 
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where,  for  simplicity,  i  and  j  are  each  assumed  tc  have  N 
types  cf  forces. 

A  2N-tuple  vector,  z  -  (x,  y)  which  satisfies  equa¬ 
tion  3.1  is  an  equilibrium  solution.  Like  many  nonlinear 
systems  cf  equations,  equations  3.1  have  more  than  one  equi¬ 
librium  point.  Geometrically,  these  equilibrium  points  are 
at  the  intersections  of  a  set  of  hypersurfaces  in  the 
2N- dimensional  space.  To  help  in  visualizing  the  geometry, 
we  can  look  at  an  example  using  a  1*1  system  as  shown  in 
figure  3.1.  In  this  case  the  hypersurf aces  simply  reduce  to 
hyperbolic  curves. 


figure  3.1  Equilibrium  Points  at  Hyperbolic  Intersections. 

The  existence  of  multiple  equilibria  makes  the  anal 
ysis  of  the  N* N  problem  very  interesting  but  difficult.  In 
chapter  IV,  some  illustrations  cn  how  the  locations  of  these 
equilibria  affect  phase  trajectories  will  be  presented. 

A  few  other  interesting  questions  arise  spontane¬ 
ously.  lor  instance,  hr  .  many  of  these  equilibria  are  there 
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Figure  3.3  Repeated  Equilibria. 


2.  Stability  and  Domains  of  Attraction 

Each  equilibrium  point  in  an  N*N  system  may  cr  may 
not  be  stable  depending  on  whether  or  not  its  equilibrium 
point  can  be  maintained.  The  property  of  neighborhood 
stability  is  important  because  it  has  a  strong  influence  on 
the  phase  trajectories.  Generally,  if  an  operating  poinc  is 
stable  (the  maximum  eigenvalue  is  negative),  then  any 
perturbation  away  from  that  point  results  in  the  system 
returning  to  the  same  point.  Conversely,  perturbations 
about  an  unstable  point  results  in  divergence  from  that 
point . 

The  notion  of  domains  of  attraction  is  also  critical 
when  determining  phase  trajectories.  Any  operating  point 
within  this  domain  cr  region  will  be  "attracted"  toward  a 
stable  eguilibrium  point.  In  short,  a  domain  of  attraction 
is  a  volume  in  the  2 N-dimen sional  space  surrounding  a  stable 
equilibrium  point.  Figure  3.4  shows  a  typical  domain  in 
which  seme  of  the  trajectories  are  shown  converging  to  a 
stable  equilibrium  point. 


//  boundary 


/i  y 


trajectory 


X  unstable  point 


stable  jjoint 


Figure  3.4  Domains  of  Attraction. 


Domains  of  attraction  are  separated  by  boundaries 
which  are  invariant  curves  in  1*1  problems  and  invariant 
hypersurfaces  in  N  *  N  problems.  A  boundary  surface  may  be 
considered  as  an  infinite  number  of  invariant  curves  placed 
side  by  side.  A  boundary  curve  is  the  locus  of  points  that 
approach  an  unstable  point  from  both  sides.  The  boundary 
line  can  be  obtained  by  backward  integration  (i.e.  using 
negative  time  in  equation  2.1)  starting  just  on  either  side 
of  an  unstable  point.  The  rationale  behind  this  method  is 
that  tc  approach  an  unstable  equilibrium,  a  point  must 
remain  exactly  on  the  boundary.  If  this  is  not  the  case, 
then  the  point  will  be  attracted  into  the  domains  and  move 


toward  a  stable  point  or  infinity.  By  performing  a  backward 
integration,  we  are  actually  retracing  the  path  taken  by  a 
point  which  previously  approached  the  unstable  equilibrium 
point.  This  method  requires  knowledge  of  the  unstable  equi¬ 
libria,  but  this  is  made  feasible  because  the  Continuation 
Methods  can  be  used  tc  find  all  equilibrium  solutions. 

B.  fINDING  THE  EQUIIIBBIOfl  SOIDTIONS 

To  obtain  a  set  cf  equilibrium  solutions,  one  has  to 
solve  eguation  3.1,  which  can  be  written  using  a  more 
compact  notation  as 

F(z)  =  0  (eqn  3.2) 

where 

F(.)  represents  the  right-hand  side  of 
equation  3.1 

z  «  (x,y) 

0  =  zero  vector 

It  is  well  known  that  numerical  techniques  for  solving 
nonlinear  equations  are  not  always  successful.  Since  equa¬ 
tion  3.2  describes  a  bilinear  system,  one  should  expect  to 
face  similiar  difficulties  when  attempting  to  solve  it 
numerically. 

Most  numerical  methods  for  root  finding  generally 
require  that  a  fairly  good  initial  guess  (z^)  be  known  so 
that  some  convergent  iteration  process 


brings  the  approximated  root  closer  and  closec  to  7  the 
desired  equilibrium  solution  or  root.  In  practice,  the 
following  difficulties  are  often  encountered  : 

(1)  The  convergence  condition  of  the  algorithm 
must  be  ensured ; 

(2)  Finding  an  initial  guess  that  is  sufficiently 
close  to  the  correct  solution  is  difficult, 
especially  for  higher  dimensions; 

(3)  Even  if  a  good  initial  guess  has  been  obtained, 
the  numerical  process  may  still  be  plagued  by 
ill-conditioning,  saddle  points,  etc.; 

(h)  Not  all  the  solutions  are  guaranteed  to  be 
found. 

1  •  Con tin ua tion  Method 

Fortunately,  the  above  problems  are  avoided  if  a 
numerical  method  called  the  Continuation  Method  [Ref.  8]  is 
used.  This  technique,  which  is  sometimes  called  The 
Imbedding  Method,  has  been  successfully  applied  in  many 
fields.  It  introduces  an  artifical  guide  which  will  channel 
the  iterates  toward  a  specific  solution.  Such  a  guiding 
principle  is  actually  a  knowledge  of  the  existence  of  a 
suitable  curve  connecting  an  initial  point  with  the  desired 
solution. 

Continuation  Method  has  significant  advantages  over 
other  numerical  techniques.  Most  importantly,  a  good  initial 
guess  is  not  necessary  and  all  the  solutions  can  be 
obtained. 

a.  Basic  Theory 

Given  the  problem  F(z)  =  0  to  solve,  the  first 
step  is  to  embed  it  into  a  homotopy  or  a  parameterized  set 
of  problems,  H  (z,t)  .  The  requirements  on  H(z,t)  are  : 


(1)  H(z,1)  =  (z)  *  0  is  the  original  problem 

(2)  H(z,0)  =  F0  (z)  =  0  has  a  trivial  or  easily 

computed  solution 

For  example,  a  homotopy  could  be  : 

H(z,t)  =  tF(z)  +  (l-t)F0(z)  ,  te[0,l]  (eqn. 3. 3) 

Using  the  above  parameterization,  the  simple 
problem  of  F0(z)  =  0  is  deformed  into  the  desired  one,  (z) 
=  0  .  This  is  done  by  calculating  the  solution  tc  the 

deformed  problem  at  each  stage  of  the  deformation.  The  exis¬ 
tence  of  a  continuous  curve  such  that  H(z(t),  t)  is  a  solu¬ 
tion  to  H(.,.)  =  0  for  all  t  £  [0,  1]  is  assumed. 

t.  Implementation 

To  actually  carry  out  the  above  continuation 
process  one  usually  differentiates  H  ( . , . )  to  form 

H(z (t) ,t)  *  0  (eqn  3.4) 

Using  eguation  3.4,  z  can  be  written  as  a  function  of  z  and 
t  as  given  in  equation  3.5.  The  function,  h(.,.)  is  prefer¬ 
ably  a  linear  function  that  can  be  integrated  numerically. 

z  =  h (z, t)  (eqn  3. 5) 

Together  with  the  initial  condition  z  (0)  =  zQ  , 

equation  3.5  is  actually  an  initial  value  problem  which  can 
be  integrated  numerically.  The  solution  at  t=1  is  then  the 
solution  to  the  original  problem  F(z)  =  0. 


2.  Algorithm  to  Obtain  2*2  Equilibrium  Problem 

A  2*2  Lanchester  problem  is  first  formulated  into  a 
Continuation  process.  It  is  followed  by  a  discussion  on  how 
the  accuracy  of  the  method  can  be  improved.  The  last  part  of 
this  subsection  includes  a  note  on  the  number  of  equilibrium 
points  in  an  N*N  problem. 

a.  Formulation 

For  the  2*2  problem,  F{z)  =0  is  explictly 


-Jl(ul*anz3*a12:4) 

+ 

"bll 23 

b  1 2  2  4 

-z2(u2*a,1z..a2.;4) 

r-> 

4» 

‘ b2 1 2  3 

b  2  2  2  4 

+ 

r  3 

" dl 1 2 1 

d  2 1 2  2 

' 2  4 ^U4+C12I1*C22Z2^ 

+ 

r4 

~ d  1 2  2 1 

d  2  2  2  2 

The  hoffotopy  is  formed  by  writing 


HjCl)  =  H2(zf0) 

+ 

1  (rl'bll23~b1224^ 

O 

II 

H  2 ( z )  =  H2(z,0) 

+ 

t(r2*b2lz3*b22z4) 

=  0 

H3(z)  =  H3(z,0) 

+ 

t(r3'dllZl'd21Z2) 

=  0 

H4(z)  =  H  4 ( z , 0 ) 

+ 

t  ^  r4 "dl 2Z1 ~  d2  2  z  2  ^ 

=  0 

where 


(eqn  5.6) 


Hj(2,0)  =  - z ^ (u^ +o ^ ^ z - ^ )  -  0 

H-,  ( z ,  0)  -  "z2^'"17  +  ‘121z3  +  ^",224^  ~  ^ 
Hj(z,0)  *  z  3  (u3+ *-  j  j z i +  ^ 2 1  2 2 ^  -  ^ 
H4(z,0)  =  _z4  ^u4  +  ci2Zl  +  C22:i2^  =  0 
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Next,  we  differentiate  equation  3.6  with  respect  to  t  and 
put  it  in  a  matrix  form 


Az  =  B 


(eqn  3.7) 


where 


Vanz3+ai:24 


allzl+bllt  a12Zl+b12t 


0 

u2+a2lVa22=4 

:n:3*dii' 

12-4^d12t 

C21Z3+d21t 

c22Z4+d22t 

VcnVc 

0 

0  VC12Z1+C22Z2 


Z  =^Z1,Z2’“3,Z4^ 


rl'bllZ3‘b12Z4 
B.  r2'b21:3'b2:-4 
r3_dllZl'd21Z2 
r4'd12Zl'd22z2 


Equation  3.7  can  now  he  integrated  numerically  using  one  of 
the  readily  available  integration  routines. 

We  have  assumed  that  the  trivial  solution  to 
H(z,  0)  =  0  has  been  previously  found. 

b.  Improving  Accuracy 

Numerical  integration  of  equation  3.7  inevitably 
produces  seme  errors  at  each  iteration.  Since  the 


[•] 


Continuation  method  relies  on  following  curves  to  arrive  at 
the  desired  solution,  it  is  esssential  that  each  iterate 
remains  close  to  the  actual  curve.  It  is  necessary  to 
include  a  way  to  correct  the  approximated  position  ty  means 
cf  a  corrector  step.  The  combination  of  integration  and 
correction  is  often  called  a  "predictor-corrector  step" 

This  process  of  prediction-correction  is  shewn 
in  Figure  3.5  where  each  integration  error  has  been  exagger¬ 
ated  for  illustrative  purposes.  The  algorithm  tc  be 
presented  later  employs  an  IMSL  routine  called  ZSCNT  for  the 
predictor  step.  Other  forms  of  curve  following  routine  can 
also  be  found  in  the  literature,  and  are  briefly  mentioned 
in  [Bef.  8]. 

c.  Trivial  Solution 

The  trivial  system  H  (z,  0)  =  0  was  chosen  to  be 

Hj(r,0)  -  -  z  ^  ( u  ^  ♦  a  ^  ^  z  ^  +  a  ^  ->  z  ^ )  -  0 

Hp(z,0)  -  ‘  (Ut+3  I  +  -)Z  , )  -  0 

L  >.  i.  1 10^4  (eqn  3>8J 

H.U,0)  *  -lslu3*cllVcnI2)  -  0 

H4(z,0)  -  -z4fu4+c12zi+C22Z2)  =  0 

In  non-degenerate  cases,  there  are  six  scluticns 
corresponding  to  eguation  3.8.  The  result  is  derived  in 
Appendix  B  which  also  deduces  the  number  of  trivial  solu¬ 
tions  for  an  N*N  problem  to  be 


i»0 


The  method  of  obtaining  the  trivial  solutions  is  given  in 
Appendix  3.  Using  a  combinatorial  identity,  N,  can  be 

l\ 

written  as 

3  1 


Figure  3.5  Integration  Path  using  Predictor-Corrector. 


Each  continuation  process  starts  from  a  trivial 
solution  z 0  and  follows  a  specific  curve  until  it  reaches 
the  equilibrium  point.  Consequently,  the  number  of  equilib¬ 
rium  points  will  also  be  Nk.  As  mentioned  in  section  IIIA, 
the  two  exceptions  are  situations  involving  inf initely-many 
and  repeated  equilibria.  Situations  involving  degeneracy 
are  discussed  in  Appendix  B. 

d.  Algorithm 

l1)  singularity  l£eatment.  In  Continuation 
Method  algorithms  [Eef.  8],  it  is  sometimes  necessary  to 
give  special  treatment  to  cases  in  which  the  curves  being 


followed  by  the  integration  routine  pass  through  a  singu¬ 
larity.  Experimentally,  it  had  been  observed  that  in  our 
problem,  the  singularity  took  on  the  form  shown  in  Figure 
3.6.  Corrective  measures  were  necessary  to  ensure  that  upon 
crossing  the  singularity,  the  large  magnitude  was  preserved 
but  the  sign  was  changed;  otherwise  the  curve  might  termi¬ 
nate  at  an  equilibrium  point  which  was  not  the  intended  one. 

In  the  algorithm,  the  presence  of  the 
singularity  is  detected  by  monitoring  the  rate  of  change  of 
the  individual  component  zi.  Once  identified,  this  fast¬ 
changing  and  large-magnitude  component  (zp)  is  monitored  at 

each  step  t  where  0  =  tn  <...<...  t.  <  t,  ,  <...<  t  .  =  1. 

u  kk+1  end 

When  zF  is  found  not  to  cross  the  singularity  and  end  up  at 
approximately  -zF  ,  the  algorithm  attempts  to  correct  this 
irregularity  by  artificially  making  zp  =  -z  before  the  next 
predictor  step  commences. 

(2)  Flowchart.  The  flowchart  for  the  algo¬ 
rithm  is  given  in  figure  3.7.  Only  the  major  steps  have  been 
shown.  The  program  listing  is  given  in  Appendix  C. 

3.  Example  and  Eesults 

Consider  as  an  example  a  2*2  problem  with  the 
following  attrition  coefficients 
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Figure  3.6  Curve  Passing  through  Singularity. 

The  trivial  solutions  are  first  computed  and  serve 
as  one  of  the  inputs  to  the  program.  The  program  obtains  the 
values  of  z  (t)  and  plots  each  component  (  (t)  j  versus  t. 

In  Figure  3.8,  the  plcts  for  t  close  to  zero  show  one  set 
of  curves  for  z^t)  starting  from  their  respective  trivial 
solutions.  The  curves  of  z^t)  versus  t  for  all  the  six 
sets  cf  eguilibrium  solutions  are  shown  in  Figure  3.9. 

A  few  interesting  features  of  the  continuation 
process  are  worth  noting.  For  example 

•  Each  trivial  solution  leads  to  different  equilib¬ 
rium  solution  and  the  integration  path  is  different 
for  each  component. 

•  All  the  curves  are  smooth  ;  one  of  the  four  curves 
may  pass  through  a  singularity.  (  see  Figure  3.9  (c) 
and  (d)  )  . 

Table  I  summarizes  the  computed  equilibrium  solu¬ 
tions.  They  are  tabulated  in  the  same  order  as  the  plots  in 


Figure  3,7  Flowchart  for 
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Figure  3.8  z  ^  (t)  Versus  t  for  Values  of  t  Close  to  Zero. 


Figure  3.9.  To  estimate  the  accuracies  of  the  results,  we 
defined  error  as 


ERROR 


+ 


+ 


-  2 


where 

Fi  =  >  2  is  t^e  computed  equilibrium 

solution  to  F(z)  =  0 

Decreasing  the  integration  step  size  in  the 
predictor  and  corrector  routines  may  reduce  the  errors  by  a 
small  amount;  but  the  increase  in  computational  effort  iray 
not  be  justifiable.  Conversely,  it  may  be  desirable  tc  cut 
down  computing  time.  Currently,  the  algorithm  performs  tne 
corrector  step  for  each  predictor  step.  If  two  or  more 


Figure  3.9  (contd.) 


TABLE  I 

Computed  Equilibria  for  X.  =  1.0,  Y  =  4.0 


Trivial  Solution 

Computed  Ecui librium  Solution,  z 

Trrcr 

(0,  0, 

0,  0) 

(1.7755,  1.2906, 

1.0033,  0.3512 

0 . 2  1 0  ~ 4 

(0,  -0. 33, 

0 ,  -0.33) 

(  34  .  3  939  ,  -04  .295  1  , 

-0.1  35  0  ,  -0 . 2775  ) 

C  .6*10'° 

(  0  ,  -0.091, 
-0.5,  0  ) 

(0.6154,  0.3345, 

3.0769,  C.9231) 

0  .  1* 10'" 

(-0.2,  0, 

0,  -1.0) 

(0.4&nC~4,  -0  .  3  9  30  , 
-16.6400,  52.9320  ) 

0  .  2  2  ••  •  1 0  "  5 

(-0.033,  0 , 
-0.3,  0) 

(-0.1357,  0.0325, 
-44.5297,  23.6147) 

0  .  5  6”10_° 

(-0.42,  0.37, 
0.25,  -0.17) 

(-45.2284,  39.1604, 

-0  .  34  97  ,  -0 . 044  85  ) 

0  .23*10"° 

predictor  steps  are  done  for  each  corrector  step,  seme 
computational  effort2  can  be  saved. 


2Saving  in  computational  effort  will  be  more  significant 
when  solving  higher  dimensional  systems. 
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IV.  PBCFERTIES  OF  THE  1*1  SYSTEM 

Tie  1*1  problem  is  the  simplest  case  in  our  model.  It  is 
nevertheless  important  for  us  to  investigate  and  understand 
its  properties.  Despite  its  relative  simplicity,  it  is  by  no 
means  uninteresting.  There  exists  many  situations  which  can 
te  realistically  and  easily  modeled  by  the  1*1  system.  For 
example,  when  the  opposing  forces  can  be  considered  as  homo¬ 
geneous,  it  is  convenient  to  use  the  1*1  model  for  analysis. 
It  is  also  useful  for  the  analyses  at  the  strategic  level 
when  the  forces  and  parameters  can  be  aggregated.  In  many 
instances,  it  seems  to  provide  insight  on  how  to  approach 
the  N*N  problem,  which  is  much  more  difficult  to  visualize. 
In  fact,  as  the  understanding  of  the  1*1  system  increases, 
there  is  a  strong  urge  to  try  to  represent  the  N*N  problem 
by  an  equivalent  1*1  problem.  The  equivalent  representation 
is  not  only  attractive  in  terms  of  its  simplicity  but  also 
its  economy  in  compu tationa 1  efforts. 

The  next  section  will  focus  on  the  relation  between 
system  asymptotes  and  stability  of  the  equilibria.  By  formu¬ 
lating  the  problem  quantitatively,  we  are  able  to  arrive  at 
some  useful  properties.  In  Section  IVB,  the  system  dynamics 
i.e.  the  changes  in  the  force  levels  are  analysed  by  consid¬ 
ering  the  phase  trajectories. 


A.  SYSTEM  ASYMPTOTES  AMD  E  QUIIIBRIUM  POINTS 


For  the  1*1  problem,  the  system  reduces  to 

x  =  -x(u  +ay)  +  r  -  by 
Y  =  -y(v  +cx)  +  s  -  dx 


(eqn  4.1) 


no 


HI'!  '  .  '  . 


An  equilibrium  condition  exists  if  r  and  s  are  chosen  such 
that  x  and  y  are  both  zero.  In  general,  there  will  be  two 
equilibrium  points  corresponding  to  two  locations  where  the 
two  hyperbolas  intersect.  The  hyperbolas  are  described  by 


r  -  bv 


u  +  ay 


s  -  dx 
v  +  cx 


(eqn  4.2) 


From  equation  4.2,  one  can  easily  deduce  the  fouc  asymptotes 
(two  vertical  and  two  horizontal)  associated  with  the  hyper¬ 
bolas.  Figure  4.1  shows  a  typical  set  of  four  asymptotes. 
They  always  cross  in  the  third  guadrant  of  the  x- y  plane  and 
do  not  depend  on  the  replenishment  coefficients.  The  rela¬ 
tive  displacements  between  the  two  horizontal  (and  also 
vertical)  asymptotes  depend  only  on  the  ratios  of  attrition 
coefficients  and  not  on  the  coefficients  themselves.  It 
turns  out  that  these  properties  of  the  system  asymptotes 
help  to  simplify  the  analysis  considerably . 


1  •  Stability  Criteria 


(xe 


/  ye) 


Considering  small  perturbations  about  an 
and  linearizing  the  equations,  we  have 


(u  +  aye)  (b  +  axe) 

Sx 

- 

(d  +  cye)  (v  +  cxe)_ 

Sy 

equilibrium 


The  characteristic  polynomial  is  simply 


Figure  4.1  System  Asymptotes. 


D(s)  =  Det  [si  -  C]  (eqn  4.3) 

where 

I  =  identity  matrix 
C  =  the  2*2  matrix  in  equation  4.3 

Hence , 

D(s)  =  [s  +  (u+aye)j£s  +  (v+cxe)  -  j\b*axe)  (d  +  cye)] 

=  s2  +  £(u+aye)  +  (v  +  cxe)  s  +  £(u+aye)  (v+cxe) 

-  (b  +  axe)  (d  +  cye)J 

Ihe  conditions  for  (xe  ,  ye)  to  be  a  stable  equilibrium, 
i.e.  for  the  roots  c£  D(s)  to  be  in  the  Left  Half  plane 
(LHP)  are  given  by 


(eqn  4.4) 


(u  +  ay  )  + (v  +  cx  )  >  0 
0  6 

(u  +  aye) (v  +  cxe)  -  (b  +  axg)  (d  +  cye)  >  0 

2.  Stability  and  Asymptotes 

five  different  ways  in  which  the  hyperbolas  can 
intercept  have  been  identified  and  their  stabilities 
accounted  for.  These  five  cases  are  shown  in  Figure  4.2  and 
each  case  will  be  elaborated  upon  subsequently. 

a.  Definitions  and  Formulation 

One  of  the  most  intriguing  facets  of  the  1*1 
problem  is  the  connection  between  the  asymptotes  and  the 
stability  cf  the  resulting  equilibria.  We  begin  the  quanti¬ 
tative  treatment  by  first  defining  the  following  ratios: 


n  -  H  n  -  d 

1  a  ’  2  "  c 


The  four  asymptotes  are  x  =  ,  x  =  -u2  ,  y  *-n  and  y  - 

-n2  •  If  we  let  the  first  equilibrium  point  be  (*el  »  Ye L) 
and  substitute  the  corresponding  r  and  s  into  the  equation 
4.2,  we  have 

Vx-Xei)  +  (xy'xelyel)  +  yl(y*yel)  =  0 

(eqn  4.5) 

n2(x'xel)  +  (xy'xelyel5  +  W2(y‘yel)  =  0 
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Next,  the  distances  between  the  asymptotes  are  dafined  as 


e  =  V  -  u 
x  2  1 


n2  -  ni 


and  f  will 


It  is  not  difficult  to  see  that 
decide  where  the  hyperbolas  intersect.  For  instance,  when  e 
>  0  and  >  0,  there  may  be  two  equilibria  in  the  first 
quadrant3  (See  case  (b)  of  Figure  4.2).  In  general,  the 
second  equilibrium  point  (Xg2,  Ye2  )  can  be  found  by  elimi¬ 
nating  y  or  x  from  equation  4.5  and  comparing  coefficients 
with  (y  -  yel)  (y  -  ye2)  and  (x  -  xgl)  (x  -  xe2)  .  The  final 
expressions  are 


x  2  =  — 
e2  e 

<yei  * 

np  - 

U1 

J 

G 

(eqn  4.6) 

v  =  _x 

e2  Ex 

<xel  * 

V  * 

nl 

For  constant  x 

el'  Xe  2 

'  yel# 

ye2,  equation 

4.6  can  be 

written 

to  represent 

two  straight 

lines  in  ex  , 

ey  plane. 

The  equations 

of  these  two  lines 

are 

e  e 

(ye2 

+  V 

y  =  x 

(xel 

+  V 

(eqn  4.7) 

e  e 

V  =  X 

*  V 

C*e2  *  Uj) 


3In  our  coi.  .ext,  the  quadrants  are  defined  by  the  asymp¬ 
totes  and  not  by  the  x,  y  axes. 
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b.  Types  of  Equilibria  and  their  Stability 

To  derive  the  different  types  of  eguilibria  and 
their  associated  stabilities,  we  make  a  transition  from  the 
i,  y  plane  into  tie  ex,  ey  plane.  Briefly,  the  basic 
approach  is  to  fix  ore  equilibrium  point  (x  .  ,  y  )  on  the 
first  quadrant  hyperbolas  and  consider  the  regions  in  the 
ex,  ty  plane  when  we  have  the  other  point  in  various  places 
of  the  x,  y  plane.  The  other  essential  step  is  to  express 
the  stability  criteria  (equation  4.4)  in  terms  of  ex,  e  , 
xel,  y  .  A  summary  of  the  results  which  are  derived 
in  Appendix  0  is  given  below  : 

(1)  When  both  equilibrium  points  are  on  the  first 
guadrant  hyperbolas  (  case  (b)  in  Figure  4.2  ),  cne 
will  be  stable  and  the  other  unstable; 

(2)  When  one  equilibrium  point  is  on  the  first 
guadrant  and  the  other  on  the  third,  both  can  be 
unstable  or  one  will  be  stable  and  the  other 
unstable  (  case  (a)  in  Figure  4.2  ); 

(3)  When  both  equilibrium  points  are  on  the  third 
guadrant  hyperbolas,  both  are  unstable  (  case  (c)  in 
Figure  4.2); 

(4)  When  there  are  infinite  number  of  equilibria  as 
in  case  (d)  in  Figure  4.2,  Gx  =  Gy  =  0  and  the  two 
sets  of  hyperbolas  merge.  Equilibria  lying  on  the 
first  quadrant  hyperbola  are  neutrally  stable  (one 
eigenvalue  equals  zero)  and  those  on  the  ether 
hyperbola  are  unstable; 

(5)  When  there  are  repeated  equilibria  as  in  case 
(e)  in  Figure  4.2,  they  are  neutrally  stable  if  the 
hyperbolas  touch  in  the  first  quadrant  ;  otherwise 
they  are  unstable. 


Most  of  the  above  results  are  embedded  within 
Figure  4.3  which  is  reproduced  from  Appendix  D  for  conven¬ 
ience.  Evidently,  hcth  the  coordinates  of  the  equilibrium 
points  (xe,  ye)  and  the  location  in  the  ex,  ey  plane  deter¬ 
mine  the  stabilities.  The  ex#  plane  has  been  subdivided 

into  a  few  regicns  each  with  distinct  stability 
characteristics. 


(•xci>  yel) 


unstable; 
^  >p7) 


e2’  >c2J 
stable 


as  above 


both 

unstable 


y„J 


el*  ' elJ 
stable; 


(-xe2’  ye2} 


unstable 
as  above 


e  =c 


<V>eP 


y  x<v3‘cl) 


Figure  4.3  The  rx,  cy  Plane. 

The  case  of  infinitely  many  equilibria  corre¬ 
sponds  to  the  origin  of  ex,  ey  plane  (  bx  =  u2  ,  n1  =  n2  ) .  The 
only  way  for  two  sets  of  hyperbolas  to  merge  is  for  their 
respective  asymptotes  to  merge.  This  case  is  a  degenerate 


instance  of  repeated  equilibria  (  case  (e)  in  Figure  4.2  ), 
which  is  shown  in  Appendix  D  to  correspond  to  operating 
points  on  the  line  ey  =  (Y  ♦  n^/CX  ♦  u  x)  as  illustrated 
in  Figure  4.3. 

As  a  corollary,  we  note  that  there  cannot  be  two 
stable  equilibrium  points  in  the  1*1  problem.  This  deduction 
can  be  made  by  referring  to  Figure  4.3.  There  is  no  region 
in  the  ex,  ey  plane  which  allows  for  this  case.  At  most, 
there  can  be  two  neutrally  stable  equilibria  which  are 
repeated.  Numerous  attempts  have  been  made  to  obtain  two 
stable  equilibria  in  the  2*2  problem,  but  in  vain.  Whether 
it  is  also  true  for  2*2  or  higher  dimensional  problems  that 
only  cne  equilibrium  may  be  stable  is  still  a  matter  of 
conjecture. 

In  Appendix  E,  the  relations  between  the  regions 
on  the  ex  ,  ey  plane  and  their  associated  stabilities  are 
verified.  Some  representative  points  on  the  ex,  ey  plane 
are  chosen  and  thedr  stabilities  checked. 


E.  SISTEH  DTNAHICS 


The  dynamics  of  a  1*1  system  are  characterised  by  its 
phase  trajectories,  which  are  curves  on  the  x-y  plane 
describing  the  history  of  the  system  as  the  time,  t, 
changes.  These  trajectories  can  be  conveniently  obtained  by 
integrating  equation  4.1  numerically. 

Needless  to  say,  being  able  to  predict  the  trajectories 
is  important,  for  it  means  that  we  know  how  our  model  of  a 
battle  progresses.  Cnee  the  factors  influencing  the  course 
of  a  battle  are  known,  appropriate  command  decisions  can  be 
introduced  to  ensure  favorable  battle  outcome.  In  Chapter  V, 
we  will  see  how  many  of  the  results  obtained  in  this  section 
can  be  used  to  rationalize  and  predict  battle  outcome. 


Some  typical  trajectories  corresponding  to  the  different 
types  of  equilibria  are  described  in  the  next  subsection. 
Besides  the  stability  which  influences  trajectories,  it  was 
briefly  mentioned  in  Chapter  III  that  domains  of  attraction 
also  affect  the  trajectories.  In  the  subsection  that 
follows,  we  will  show  specific  examples  of  the  way  to  deter¬ 
mine  the  domains  by  finding  their  exact  boundaries. 

1.  Ira  iectories 

Two  methods  of  establishing  the  trajectories  from  a 
given  initial  condition  will  be  described.  The  brute-force 
method  which  has  been  mentioned  uses  numerical  integration. 
The  ether  method  which  often  provides  better  insight,  is 
more  graphical.  The  graphical  method  is  based  on  a  few  very 
simple  rules  to  predict  the  gross  behavior  of  a  trajectory. 
Some  cf  these  rules  are  listed  below  : 

(1)  A  stable  point  "attracts";  unstable  point 

"repels" ; 

(2)  Points  on  either  side  of  a  boundary  move  into 

their  respective  domains; 

(3)  For  large  (x,  y)  ,  trajectories  are  governed  by 

the  Lanchester  "linear  law"; 

(i»)  Points  near  the  hyperbolas  can  be  easily 

•  • 

analyzed  by  noting  the  signs  of  x  and  y. 

As  an  example  of  using  the  graphical  method  to 
determine  trajectories,  consider  a  region  around  an  unstable 
equilibrium  point  on  the  first  quadrant  hyperbola.  The  whole 
picture  of  the  phase  trajectories  (sometimes  called  phase- 
plane  portrait  [Ref.  9]  )  can  be  put  together  in  a  logical 
fashion  by  using  those  simple  rules.  Since  this  equilibrium 
point  is  unstable,  trajectories  will  be  expected  to  diverge 
from  it.  As  an  unstable  equilibrium  point,  it  will  have  a 
boundary  line  passing  through  it.  Initial  conditions  start 


from  each  side  give  rise  to  different  trajectories.  Next,  we 
determine  the  signs  cf  x  ,  y  on  both  sides  of  each  hyperbola 
as  indicated  in  Figure  4.4  where  only  one  intersection  is 
shown. 


X<0 

o 


^hyperbola 


Figure  4.4  Analytical  method  of  predicting  trajectories. 


Note  how  predictable  these  trajectories  are.  If, 
for  some  reasons,  the  exact  trajectories  are  required,  we 
can  resort  to  the  brute-force  method.  The  methods  are  obvi¬ 
ously  complementary  in  nature.  The  advantages  of  the  brute- 
force  method  are  accuracy  and  simplicity.  In  Figure  4.5,  a 
typical  computer  plot  consisting  of  ten  trajectories  is 
shown.  The  program  which  produces  the  plot  is  included  in 
Appendix  F. 

Beferring  to  Figure  4.5,  the  trajectories  cross  the 
hyperbolas  and  move  asymptotically  along  a  common  curve 


figure  4.5  Computer  Plot  of  Trajectories. 

lying  between  the  hyperbolas.  This  same  property  is  exhib¬ 
ited  by  other  cases.  Even  the  special  case  with  no  hyper¬ 
bolic  intersection  has  been  found  to  behave  similarly  as  can 
seen  in  Figure  4.6. 

Our  ability  to  determine  the  trajectories  and 
present  them  vividly  is  partly  due  to  fact  that  two- 
dimensional  pictures  can  be  easily  drawn  and  visualized.  For 
dimensions  higher  than  the  third,  it  is  impossible  to  visu¬ 
alize  trajectories;  however,  the  notion  of  trajectories  can 
be  conceptually  extended  to  n-dimensional  space.  Thus,  it 
seems  likely  that  in  the  higher  dimensional  systems,  trajec¬ 
tories  cross  hyperscrfaces  and  move  along  a  common  asymp¬ 
totic  curve  analogous  to  that  in  the  1*1  system.  Further 
studies  are  required  before  this  behavior  can  be  confirmed. 


y 


Figure  4.6  Trajectories  when  Hyperbolas  do  not  Intersect. 

2.  Boundaries  of  Domains  of  Attraction 

In  Chapter  III,  the  idea  of  the  domains  of  attrac¬ 
tion  was  briefly  discussed.  In  an  n-dimensional  space,  such 
a  domain  is  a  region  or  volume  in  which  all  initial  points 
come  under  similiar  influence.  When  domains  exist,  there 
will  be  boundary  surfaces  which  can  be  thought  cf  as 
collections  of  invariant  curves  passing  through  unstable 
equilibria. 

For  a  1*1  problem,  domains  and  boundaries  are  net  at 
all  abstract.  In  the  last  subsection,  they  have  been  shown 
to  affect  trajectories.  Recall  that  in  Chapter  III,  we 
mentioned  a  simple  and  yet  effective  way  of  finding  the 
boundary  curves  and  establishing  the  domains  in  the  x-y 
plane.  Examples  on  tie  use  of  backward  integration  to  obtain 
boundary  curves  are  now  presented. 
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a.  Boundary  Curve  through  an  Unstable  Point 

Starting  from  an  unstable  point,  we  apply  small 
perturbations  in  both  directions  perpendicular  to  an  eigen¬ 
vector  associated  with  the  most  positive  eigenvalue  and 
integrate  backward  in  time  (in  the  computer  program,  this  is 
easily  done  by  employing  negative  time  steps  for  integra¬ 
tion)  .  The  result  is  a  smooth,  invariant  curve  which  is 
exactly  the  boundary  or  the  so-called  separatrix  like  the 
one  shown  in  Figure  h.7. 


Boundary 


Figure  h.7  Boundary  Curve  through  an  Unstable  Point. 

To  verify  that  the  curve  is  indeed  the  boundary, 
two  initial  points  are  chosen  just  off  the  curve  (e.g  A,  B 
in  Figure  4.7) .  if  we  forward  integrate  from  these  two 
points,  they  move  into  different  domains  as  indicated  in  the 
same  diagram.  Appendix  G  contains  a  Fortran  program  that 
does  the  backward  integration  and  plots  the  boundary  curve. 


3.  Eounda ry  Curve  between  Two  Hy perbo las 

Boundary  curves  do  not  necessarily  pass  through 
unstable  points.  Backward  integration  methods  can  also  be 
used  if  a  boundary  exists  but  there  is  no  unstable  equilib¬ 
rium  point  to  serve  as  the  starting  point  of  integration. 
This  is  best  illustrated  by  considering  the  case  of  both 
equilibria  cn  the  first  quadrant  hyperbolas.  In  this  case, 
there  is  no  equilibrium  point  in  the  third  quadrant;  never¬ 
theless  a  boundary  does  exist  between  the  third-quadrant 
hyperbolas.  The  existence  of  the  boundary  is  visible  by 
simply  ccnsidering  the  signs  of  x  and  y  on  both  sides  of  the 
hyperbolas.  In  figure  4.8,  the  signs  of  x  and  y  and  also  the 
directions  of  some  typical  trajectories  are  depicted. 


II  • 
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Figure  4,8  Existence  of  Boundary  Between  Two  Hyperbolas. 

To  obtain  the  exact  boundary,  choose  a  point  clcse 
to  a  hyperbola  and  cn  lower  part  of  the  hyperbolas  (e.g. 
point  P  or  Q  in  Figure  4.8)  and  integrate  backward.  The 
result  is  a  boundary  curve  as  shown  in  Figure  4.9. 


hyperboles 


Figure  4.9  Exact  Boundary  Curve  Between  Two  Hyperbolas. 


4 .  Summary  of  t  he  1*1  Problem 


We  have  seen  the  close  relation  between  system 
asymptotes  and  stabilities.  Through  the  use  of  newly  defined 
variables  e  and  e  ,  the  stability  of  different  types  of 

x  y 

equilibria  has  been  derived.  Five  cases  have  been  identi¬ 
fied,  and  they  correspond  to  the  types  of  intersections  on 
the  x-y  plane.  For  example,  if  both  the  equilibria  are  found 
on  the  third  quadrant  hyperbolas,  then  we  know  that  they 
will  be  unstable. 

Two  methods  cf  establishing  the  trajectories  have 
been  described  in  this  chapter.  These  two  methods  complement 
each  other  and  the  choice  depends  on  our  requirements.  The 
dynamics  of  the  system  are  characterized  by  the  trajecto¬ 
ries,  which  as  we  have  seen  are  very  predictable.  These 
trajectories  are  influenced  by  the  stabilities  of  equilibria 
and  domains  of  attraction  which  are  separated  by  boundary 
curves.  A  simple  way  of  plotting  the  boundary  curves  has 
also  been  presented  along  with  specific  examples. 


The  results  derived  in  this  chapter  will  be  applied 
in  the  next  chapter.  The  knowledge  of  the  system  dynamics 
and  how  they  are  affected  by  stability  and  other  parameters 
will  enable  us  to  analyze  changes  in  force  levels  as  the 
tattle  progresses. 


V.  STRATEGY  FOB  INITIAL  FORCE  COMMITMENT 

In  the  last  two  chapters,  emphasis  has  been  placed  on 
establishing  the  mathematical  framework  of  the  system 
dynamics  and  stability.  In  this  chapter,  we  examine  some 
model  operational  problems  that  are  related  to  stability  and 
dynamic  considerations. 

One  of  the  major  command  decisions  that  has  to  be  made 
during  a  build-up  period  of  a  war  pertains  to  initial  force 
commitment.  A  good  strategy  calls  for  a  balance  between 
initial  deployment  and  reserves.  In  practice,  a  multitude  of 
factors  have  to  be  considered  before  deciding  on  a  partic¬ 
ular  commitment.  The  approach  in  this  chapter  provides  us 
with  a  set  of  mixed  strategies  but  does  not  consider  intan- 
giable  factors  like  world  politics,  national  economy, 
survival  factor  and  so  on. 

Stability  has  been  shown  to  effect  trajectories  which  in 
turn  effect  battle  outcome.  Recall  from  Chapter  IV  that 
there  are  some  trajectories  which  represent  speedy  and 
complete  annihilation  of  one  force;  hence  it  seems  reason¬ 
able  that  the  side  that  is  tipped  to  win  the  battle  will 
want  to  operate  on  an  unstable  trajectory.  But  to  what 
extent  can  one  exploit  the  stability  behavior  of  the  system 
to  influence  battle  outcome?  Obviously  there  will  be  prac¬ 
tical  limitations;  an  important  one  of  these  is  total  avail¬ 
able  resources. 

A.  EBOBLEM  STATEMENT  AND  APPROACH 

The  problem  statement  is  as  follows  : 

Given  total  defense  resources  Qx  ,  Qy  for  x  and  y 
respectively,  what  is  the  optimum  set  of  strategies  for 
initial  force  ccnmitment,  X  and  Y? 


W €  begin  by  treating  this  as  a  1*1  problem  at  the  stra¬ 
tegic  level.  The  dynamics  of  the  problem  are  thus  governed 
by  eguation  4.1.  Both  sides  are  assumed  to  operate  initially 
at  equilibrium  with  constant  replenishment  rates  given  by 


r  =  X(u  +  aY)  +  bY 
s  =  Y(v  +  cX)  ♦  dX 


(eqn  5.1) 


Since  both  sides  have  limited 
the  replenishment  rates  versus 

Figure  5.1,  where  Q  =  rT  and  Q 

x  x  y 


defense  resources  Q  ,  Q  , 
time  may  be  as  shown  in 
-  sTy. 


replenishment 


t 

Figure  5. 1  Beplenishaent  Versus  Time. 

The  next  step  is  to  select  some  suitable  form  of  payoff 
function  which  is  to  be  optimixed  for  a  certain  choice  of  X 
and  Y.  The  payoff  function  (from  X  to  Y)  has  been  chosen  to 
be 

A (X , Y)  =  Ly  -  Lx 

where  Lx,  Ly  *  Total  losses  for  x,  y  at  battle  termination 
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As  each  side  runs  out  of  resources  at  different  times, 
the  simulation  is  conducted  in  stages.  The  total  losses  are 
determined  by  simulating  the  dynamics  of  the  system  until 
one  of  the  force  levels  drops  to  ten  percent  of  its  total 
resources,  Q. 

If  X  and  y  are  assumed  to  he  chosen  from  a  finite  set  of 
values,  then  for  each  pair  (X,Y),  one  A(X,Y)  can  be 
obtained.  A  payoff  matrix  can  be  formed  and  the  problem  can 
be  treated  as  a  two-person  game.  Based  on  the  mirimax 
theorem,  there  exists  a  set  of  optimal  mixed-strategies  and 
one  convenient  way  of  finding  them  is  through  the  use  of 
linear  Programming. 

It  is  perhaps  worth-noting  that  the  approach  is 
computation-  oriented.  It  has  been  made  feasible  by  the 
availability  of  high-speed  computers  and  efficient  software 
for  numerical  computations. 

E.  MULTISTAGE  BATTLE 

Using  the  above  approach,  the  entire  battle  can  be 
divided  into  three  stages,  namely 

(1)  Both  r  and  s  are  nonzero 

(2)  One  of  the  r  or  s  eguals  zero 

(3)  Both  r  and  s  are  zero 

1  •  Stage 

This  stage  will  be  the  period  from  outbreak  of  war 
to  the  time  (T^  when  one  side  runs  out  of  resources.  It  is 
also  possible  that  x  <0.  1Q^  or  y  <0. 1Qy  before  T1  is 
reached,  in  which  case  the  battle  is  over.  In  general,  this 
period  T.  can  be  written  mathematically  as 


During  this  stage,  the  dynamics  of  the  system  is 
given  by  the  familiar  1*1  system 


x  =  -x(u  +  ay)  -  by  +  r 
y  =  -y(v  +  cx)  -  dx  +  s 


(eqn  5.2) 


Khen  this  1*|  system  is  integrated,  just  as  in 
Chapter  IV,  the  resulting  trajectories  behave  similiarly. 
However,  there  is  a  major  difference.  Now,  we  no  longer 
have  unlimited  defense  resources,  and  this  stage  will  not 
last  forever.  It  implies  that,  unless  Qx  or  Qy  is  extremely 
large4  ,  trajectories  which  reflect  quick  annihilation  of 
one  of  the  forces  are  rare.  In  general,  Tx  and  Ty  are  given 
by 


T 


x 


-  X 


r 


s 


If  one  of  the  force  levels  drops  to  less  than  ten 
percent  of  Qx  or  Qy ,  the  battle  is  arbitrarily  considered 
over  and  the  losses  are  calculated  as  in  Figure  5.2.  Ihc 
finish  time  (FIKTIM)  is  simply  t,  the  time  when  x  <  0.  lQx  or 
y  <  0.lQy. 

2  •  Stage  2 

Since  either  x  or  y  can  run  out  of  reserves  first, 
the  dynamics  of  stage  2  are  governed  by  either  equation  5.? 
cr  5.4  respectively. 


*QX  or  Q  may  be  very  large  if  x  or  y  is  backed  by  a 
superpower  vh6  is  fglly  committed  to  provide  military  aid. 


Figure  5.2  losses  at  Stage  1. 


x  *  -x(u  +  ay)  -  by 
y  *  -y(v  ♦  cx)  -  dx  +  s 


(eqn  5.3) 


x  =  -x(u  +  ay)  -  by  +  r 
y  =  -y(v  +  cx)  -  dx 


(eqn  5.4) 


Unless  the  battle  ends  earlier,  this  period  will 
last  for  T 2  which  is  given  by 

Tz  =  Max  |TX>  Ty|  -  Tj 

During  this  period,  the  trajectory  will  be  different  from 
that  in  stage  1.  This  is  because  when  r  =  0  or  s  =  0,  cr.e  of 
the  hyperbolas  is  shifted  so  as  to  cross  the  origin  and  we 
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have  different  equilibrium  points.  The  trajectory  will  now 
be  influenced  by  the  new  equilibrium  point.  This  is  illus¬ 
trated  in  Figure  5.3. 

Calculations  of  the  losses  are  more  complex  than  in 
stage  1  since  there  are  now  two  cases  to  deal  with  i.e.  r  = 
0  or  s  =  0.  The  procedure  is  shown  in  Figure  5.4. 

3 .  Stage  3 

If  the  battle  enters  stage  3  without  either  x  < 
0.1QX  or  y  <  O.IQy  then  the  dynamics  will  be  dominated  by 
attritions  since  r  =  s  =  0.  Equation  5.5  is  now  used  for 
integration . 


x  =  -x(u  +  ay)  -  by 
y  =  -y(v  +  cx)  -  dx 


(eqn  5.5) 


Again,  the  trajectory  will  have  to  change  because 
now  both  hyperbolas  pass  through  the  origin.  This  is  illus¬ 
trated  in  Figure  5.5  where  we  show  how  the  intersection  at 
stage  2  has  changed,  losses  and  FINTIM  are  calculated  in 
accordance  with  the  procedure  in  Figure  5.6. 


C.  MIXES  STRATEGIES 


The  range  0  to  Qs  for  both  X  and  Y  can  be  subdivided 
into  m  force  levels.  There  are  m*m  pairs  of  X  and  Y  and 
corresponding  number  of  payoffs,  A(X,Y).  We  thus  have  an  m*m 
payoff  matrix  having  elements  A  (X,  Y)  -  Figure  5.7  gives  a 
pictorial  representation  of  this  two-person  game. 


5In  the  actual  program,  one  may  wish  to  restrict  the 
range  of  X  and  Y  to  interval  (0.2Q  -  0.75Q)  to  reflect  prac¬ 
tical  limitations  in  initial  force  deployment. 


gure  5.3  Trajectories  and  Hyperbolic  Intersection  I^n 


Is  r  =  0 
wo  r  s  =  0  ? 


Integrate 
I  eqn  S . 3 


x<0 . IQ  or 
s.  v<  0  .  IQ  ,? 

x  y  / 


Integrate 
eqn  5.4 


x<  0 . IQ  or 
x 

.  y<  0 . 1 Q  ? 

X  v/ 


Lx  “  Qx  *  x 

Ly  *  s(Tx  +  t)  +  Y 

A (X , Y)  *"  Ly  -  Lx 

FINTIM  =  t  +  T 


I-x  =  r(T  +  t)  +  X 
L  =  Q  -  y 
A ( X , Y )  =  Ly  -  Lx 
FINTIM  =  t  +  T 

v 


Figure  5.4  Losses  at  Stage  2. 

In  the  last  section,  the  procedure  for  computing  A  (X,Y) 
has  been  described.  A  simple  program  can  be  written  to 
compute  each  element  of  the  payoff  matrix.  One  such  program 
is  given  in  Appendix  B. 
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Trajectories  and  Hyperbolic  Intersect  ion  During  Stage 


Figure  5.6  losses  at  Stage  3. 

There  are  a  few  ways  of  presenting  and  interpreting  the 
payoff  matrix.  A  normal  practice  is  to  present  it  in  tabular 
form  and  consider  only  pure  strategy.  Alternatively,  a  plot 
of  A(X,Y)  as  a  function  of  X  and  Y  could  be  obtained.  When 
using  pure  strategies,  it  has  been  observed  that  the  game 
does  not  always  have  a  saddle  point  [Ref.  10]  and  it  would 
be  better  to  use  mixed  strategies. 

In  mixed  strategies,  x  and  y  may  play  all  their  strat¬ 
egies  in  accordance  with  a  certain  set  of  pro babilities. 
Although  in  our  situation,  x  and  y  can  only  play  once,  the 
same  concept  of  mixed  strategies  is  still  useful.  If  we  let 
p.  and  q,  be  the  probabilities  by  which  x  and  y  select  their 
ith  and  jth  pure  strategies  respectively,  then 
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Figure  5.7  Payoff  Matrix. 


Epi  ■  Bi  ■ 1 

i  j 


In  addition  the  (i,j)th  entry  of  the  payoff  matrix  be 
denoted  by  a^j,  the  probabilities  can  be  represented  by  the 
matrix  below 
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Ihe  optimal  mixed  strategy  is  based  on  the  tnir.imax 
criterion.  Mathematically,  x  and  y  select  p.  and  g.  which 


will  yield  U  and  V  as  given  equation  5.6  and  equation  5.7 
respectively. 
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(eqn  5.6) 


(eqn  5.7) 


Appendix  I  describes  how  the  problem  of  solving  for  the 
optimal  values  of  pi  and  g  can  be  put  into  linear  program¬ 
ming  form.  The  program  given  in  Appendix  H  also  computes 
this  optimum  set  of  solution  in  addition  to  obtaining  the 
payoff  matrix. 

The  concept  of  mixed  strategies  is  quite  intuitive  if  a 
game  is  to  be  played  repeatedly.  But  since  we  are  using  it 
to  provide  us  with  an  optimum  set  of  probabilities  of 
selecting  the  pure  strategies,  some  interpretation  is 
required.  Although  the  optimum  mixed  strategies  have  been 
obtained,  a  pure  strategy  still  has  to  be  selected  and  used. 
However,  it  is  important  that  the  selection  process  should 
be  random6  according  to  the  optimized  probabilities 
obtained. 

One  simple  but  valid  statistical  procedure  [Ref.  11]  to 
select  a  pure  strategy  from  a  set  of  mixed  strategies  is  to 
first  plot  the  probability  distribution  function.  A  random 
number  generator  is  then  used  to  generate  a  number  between 
zero  ard  one.  The  corresponding  value  of  the  strategy  could 
then  be  selected.  This  procedure  is  shown  in  Figure  5.8. 


6Tbe  selection  process  must  be  random  otherwise  the 
opponent  can  select  a  strategy  to  improve  his  outcome. 


P(  X  <  x  ) 


a* 
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rondom 

number 
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strotegy 
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Figure  5.8  Obtaining  Pure  Strategy  from  Mixed  Strategies. 

0.  EXAMPLE  OS IMG  KOREAN  BAB  DATA 

One  cf  the  main  objectives  of  using  actual  historical 
data  in  a  model  is  for  validation.  It  is  important  that  the 
results  obtained  using  the  model  should  at  least  be  consis¬ 
tent  with  actual  events.  The  Korean  Mar  has  been  chosen 
because  there  was  a  clear-cut  victor  during  the  initial 
phase  of  the  war.  He  consider  the  period  when  only  North 
Korea  and  Republic  of  Korea  (South  Korea)  were  involved. 

Before  the  entire  simulation  can  be  carried  out,  the 
actual  force  strategies,  fighting  ability,  weapon  state, 
etc,  have  to  be  transformed  into  familiar  quantities  and 
parameters  such  as  Qx,  Qy,  X,  Y,  a,  b,  c,  d,...,  and  sc  cn. 
This  transformation,  together  with  some  background  data  on 
the  Korean  Kar  are  given  in  Appendix  J. 

1 .  Results  and  Discussions 

First  we  examine  the  resultant  trajectories  during 
the  three  stages  of  the  battle  which  are  shown  in  Figure 
5.9.  The  simulation  uses  the  X  and  Y  which  correspond  to 
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the  actual  initial  deployment  ty  both  North  and  South  Korea 
respectively.  Clearly  we  see  that  the  victor  is  x,  as  it 
was  in  history.  The  result  of  the  simulation  also  shows  the 
three  stages  explained  in  the  last  section.  Note  that  the 
trajectories  for  the  first  and  second  stages  are  curtailed 
because  both  sides  run  out  of  war  reserves.  The  implication 
is  that  in  practice,  the  kind  of  trajectories  leading  to 
large  and  rapid  changes  in  force  levels  are  rather  rare. 

However,  the  effect  of  instability  on  battle  outcome 
is  borne  out  by  experimenting  with  the  directions  of  pertur¬ 
bations.  Consider  the  case  in  which  x  (North  Korea)  fixes 
the  initial  force  and  y  (South  Korea)  varying  the  initial 
force  levels  around  the  equilibrium  point.  In  Figure  5.9, 
these  perturbed  points  are  denoted  by  points  A  to  D  spanning 
across  the  boundary  separating  the  domains  of  attraction. 
From  our  understanding  of  the  stability  and  system  dynamics 
each  perturbation  will  give  rise  to  different  trajectory  and 
payoff  at  the  end  of  the  simulation.  Clearly,  y  will  want  to 
operate  at  the  perturbed  points  A  or  B  rather  C  or  D  since 
the  former  will  result  in  the  trajectory  for  stage  one  to  be 
in  a  decreasing  x  direction.  Table  II  shows  the  variation  in 
the  payoff  as  the  perturbation  point  changes.  When  the 
perturbed  points  are  at  A  or  E,  the  payoffs  to  x  are  less 
then  those  for  points  C  or  D.  Thus  we  have  seen  how  an 
unstable  system  can  be  used  to  inflict  heavier  losses  on  the 
opponent.  The  more  unstable  a  system  gets,  the  more  signif¬ 
icant  will  be  the  effects  of  initial  perturbation  which  are 
manifested  by  initial  victory  and  element  of  surprise.  Since 
some  systems  with  large  aimed-fire  coefficients  tend  to  be 
highly  unstable,  we  can  expect  this  effect  to  be  most 
pronounced  in  battles  involving  high-technology  and  highly- 
lethal  weapons. 

The  payoff  matrix  and  optimal  probabilities  p*  and 
q*  are  shown  in  Table  III.  The  results  suggest  that  the 


TABLE  II 

Effect  of  Different  Perturbations  on  the  Payoff 


Location  in  Co-ordinates  of  Payoff  to  X 
Figure  5.9  Perturbed  Point 

A  (6.7,  5.05)  2.39 
B  (6.7,  3.025)  2.41 
C  (6.7,  2.975)  2.45 
D  (6.7,  2.95)  2.47 


figure  5.9  Trajectory  for  Korean  War. 
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North  Koreans  should  use  large  initial  deployment.  In  the 
actual  war.  North  KoLea  actually  deployed  almost  all  of  its 
regular  force  and  within  a  few  Jays  captured  Seoul,  the 
capital  cf  South  Korea.  The  payoff  matrix  also  shows  that 
no  matter  which  strategy  is  chosen  by  South  Korea,  it  is 
bound  to  suffer  much  more  losses  than  North  Korea.  Again, 
this  is  in  agreement  with  history  since  it  is  an  accepted 
fact  that  without  US  intervention,  there  would  be  r.o  South 
Korea  today. 

So  far  in  the  example,  we  have  always  con  siuered  the 
situation  in  which  the  equilibrium  point  (X,Y)  determines 
the  replenishment  rates  as  given  in  equation  5.1.  It  is 
interesting  to  investigate  the  effect  on  the  payoff  when  the 
initial  operating  point  is  at  some  other  location  cthtr  than 
(X,Y).  let  the  new  initial  point  be  at  (K1,Y1)  and  consider 
a  case  where  X^  is  kept  equal  to  X  ar.d  only  ] i  is  varied. 
(X,Y)  has  been  cLosen  to  be  (6.7,3. 0).  In  Figure  5.10,  three 
trajectories  corresponding  to  Y^  at  2.5,  4.0,  5.0  are  shown 
together  with  the  hyperbolic  intersection  during  stage  one. 
Basically,  the  trajectories  correspond  to  the  three  stages 
of  simulation  as  before  and  x  is  still  the  victor.  However, 
fcoth  the  payoff  (Ly-Lx)  and  finish  time  ace  slightly 
different  from  operating  at  (X,Y).  Table  IV  stows  tint  y 
inflicts  mere  losses  on  x  when  operating  at  Y1  above  the 
boundary  curve  rather  than  at  (X,Y),  but  in  doing  so  y  is 
defeated  faster.  Thus  depending  on  his  mission,  a  commander 
can  choose  to  lengthen  the  battle  or  inflict  more  casualties 
on  his  opponent  by  choosing  a  suitable  operating  point  which 
nay  be  other  than  an  equilibrium  point. 
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Figure  5.10  Initial  Operating  Points  at  Non-egailibrium  Points 


TABLE  IV 

Effect  of  Operating  at  Non-equilibrium  Points 


Y1 

(Ly-Lx) 

Finish  Time 
(FINTIM) 

Remarks 

2.5 

2.475 

0.323 

Below  boundary 

3.0 

2.435 

0.313 

At  equilibrium 

4.0 

2.353 

0.293 

Above  boundary 

5.0 

2.305 

0.283 

Above  boundary 

VI.  COHCIPSIONS  AND  RECOMMENDATIONS 

A.  CONCIUSIONS 

This  thesis  has  covered  a  number  of  subjects  which  are 
based  on  the  generalized  Lanchester  Model.  The  first  part  of 
the  results  has  to  do  with  finding  the  equilibrium  points  in 
the  N*N  system.  The  Continuation  Methods  have  been  found  to 
be  suitable  for  this  purpose.  The  advantages  of  the 
Continuation  Methods  over  numerical  techniques  are  numerous 
and  important  to  our  understanding  of  the  non-linear  set  of 
equations.  The  method  finds  all  the  equilibrium  solutions 
accurately  and  does  not  need  good  initial  guesses.  An 
example  to  compute  the  equilibrium  solutions  of  a  2*2  system 
is  presented  along  with  a  way  to  treat  singularity  problem. 

The  derivations  and  interpretations  of  the  relations 
between  stability  and  system  parameters  form  the  next  major 
portion  of  the  thesis.  By  considering  the  simpler  1*1 
problem,  a  few  interesting  conclusions  have  been  reached, 
namely 

(1)  Eoth  the  system  asymptotes  and  equilibrium  points 
are  intrinsic  to  a  system  in  equilibrium.  The  locations 
cf  the  equilibrium  points  on  the  x-y  and  cx,  ey  planes 

completely  characterize  their  stabilities;  e  and  e 

x  y 

are  the  differences  in  the  system  asymptotes; 

(2)  The  dynamics  of  a  system  are  characterized  by  the 
phase  trajectories  which  represent  the  ways  a  battle 
progresses.  Besides  stability,  the  domains  of  attrac¬ 
tion  also  influence  the  trajectories.  The  -oundary 
curves  which  separate  these  domains  can  be  ascertained 
by  graphical  or  backward  integration. 


The  last  portion  of  the  thesis  integrates  the  concept  of 
equilibrium  stability  and  system  dynamics.  It  relates  these 
theoretical  concept  to  operational  problems.  Two  operational 
issues  are  addressed  namely,  (1)  the  effect  of  varying  X  and 
Y,  the  initial  force  deployment  on  battle  outcome,  (2)  the 
exploitation  of  stability  to  influence  battle  outcome.  A 
methodology  which  combines  multistage  battle  simulation  with 
two-person  game  has  been  employed  and  the  conclusions  are 

(1)  Initial  force  deployment,  X  and  Y  can  be  optimized 
by  finding  a  set  of  mixed  strategies.  A  suitable  pure 
strategy  can  then  be  selected  from  the  mixed  strat¬ 
egies; 

(2)  Instability  can  and  should  be  used  to  shape  the 
course  of  battle  and  its  outcome.  This  is  particularly 
true  in  highly  unstable  warfare  which  is  normally  asso¬ 
ciated  with  large  aimed-fire  attritions.  Unless  defense 
resources  are  extremely  large,  it  is  not  possible  to 
completely  reverse  the  outcome  of  a  lopsided-battle 
where  one  side  is  much  stronger  than  the  other. 

As  far  as  military  commanders  are  concerned,  the  above 
conclusions  suggest  two  things.  Firstly,  depending  cn  the 
relative  strengths,  it  is  not  necessarily  true  that 
deploying  the  largest  possiLle  force  will  bring  victory, 
reduce  loss  or  even  buy  time.  There  is  an  optimum  way  of 
deploying  available  forces.  Secondly,  if  a  war  involves 
large  aimed-fire  attritions  due  to  weapons  like  aircraft, 
missiles,  tanks,  artillery,  naval  bombardment,  etc.,  then 
initial  victory  which  could  perhaps  be  achieved  through  a 
preemptive  strike  certainly  affects  battle  outcome. 


£.  RECOHHEHDATICHS 


1  •  Transforma  tior  of  N  *N  Problem  into  the  1_*  ±  Problem 

It  has  been  mentioned  that  the  simplicity  of  the  1*1 
problem  can  be  attributed  to  the  simpler  mathematics 
involved  and  our  ability  to  draw  and  visualize  two- 
dimensional  pictures.  Despite  its  simplicity,  it  does  share 
many  of  the  properties  with  the  N*N  problem.  Considering 
these  factors,  it  seems  logical  that  an  attempt  should  be 
made  to  find  the  1*1  equivalent  to  the  N*N  system.  Another 


reason 

is  that  there  is 

much  to 

be 

gained 

in  terms 

of 

savings 

in  computational 

ef  fort 

by 

going 

to  the 

1*1 

equivalent. 

Of  course,  the  " 

equivalent" 

system 

will  not 

be 

expected  to  be  identical  to  the  N*N  problem  in  every  aspect. 
One  can  only  hope  that  it  is  equivalent  in  some  sense,  for 
example 


(a)  Preservation  of  stability  characteristics  and 
dynamics; 

(b)  Preservation  of  mixed  strategies. 

One  way  of  transforming  the  N*N  system  into  the 
equivalent  1*1  system  is  t o  equate  losses  in  Loth  systems. 
The  equivalent  system  parameters  are  obtained  by  using  rela¬ 
tions  such  as 


where 


aX 


Y 

eq  eq 


EE'hjVj 

i  j 


XV  r 

i 9  =  coordinates  of  the  equilibrium  point  in 

the  N*N  problem 
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She  results  of  the  preliminary  studies  suggest  that 
this  method  of  transformation  can  preserve  some  stability 
characteristics.  The  possibility  of  using  the  equivalent 
system  to  obtain  the  mixed  strategies  should  not  be 
dismissed  until  further  studies  have  been  conducted. 

2.  Time  Variable  Re  pie  n ish  ment  Coef f icients 

The  replenishment  rates  used  in  the  thesis  have 
always  been  assumed  to  be  constant.  In  actual  wars,  constant 
replenishment  rates  may  not  be  used  by  either  side;  at 
times,  it  may  not  even  be  possible  to  do  so.  It  would  there¬ 
fore  be  interesting  to  study  the  cases  which  involve  time- 
varying  replenishment  rates  r(t).  The  choice  of  r(t),  for 
example  periodic,  non-periodic,  ramp,  etc.  will  depend  on 
how  well  it  represents  practical  replenishment  rates. 
Whether  a  mathematical  tool  can  be  found  to  cope  with  the 
added  complexity  also  has  to  be  considered. 


C  APPENCIX  A 

C  SOLVING  FOR  OPTIMUM  VECTORS  AND 

C 

c 

c 

£****44444 44*44  <44*4**4***44*444***4*444**44**444**4*4*4 


C 

C 

c 

c 

c 

c 

c 


THIS  PROGRAM  COMPUTES  THE  OPTIMUM  FORCE 
OISTR  JEUTICN,  X*  GIVEN  X,Y,AND  THE  ATTRITION 
CCEFF  IC  IENTS. PARAMET ERS  K1.K2  AND  CONFLICT 
MATRIX  ARE  ALSO  QETAINEi). 

THIS  IS  AN  INTERACTIVE  PROGRAM. EEF ORE  EXECUTION 
i CHECK  TFE  VALUES  OF  ATTRITION  COEFFICIENT 
MATRICES  A  A , B  6  ,  CC  ,  C  AND  t Ui , U 2 , V 1  ,  V2 J 
IN  THE  FILE  L  PD  FORTRAN. THEN  ENSURE  YCU  HAVE 
ANOTHER  FILE  CALLED  CCLIB  EXEC. TO  EXECUTE  THE 
PROGRAM, ENTER  "CO LIB  LPD"  ANU  FOLLCw  INSTRUCTIONS. 
ON  INPUT  OF  X  AND  Y  HILL  CREATE  THE  A  MATRIX 
AND  RUN  THE  LINEAR  PROGRAM  ZX4LP. 

£**4**44444444  4 4 4 4* <**4  4*4  4 4* 44 *******  4***4 44* 4 4*** 4** *4 

C 

A(Ml+M2*2,N*Nl+2 ) 

E { Ml +M2  i 
UN) 

PSOL (N) 

OSOL  (M1  +  M2J 

RWII  A*(Ml*M2  +  4)+2*(N*MlJ  1 
IW<2*(N*M1J  ♦  3J 

REAL  At  4,  <j,B(2)  ,C(3  J ,  S  ,  PSOL  (3)  ,  CSOL  (  2) ,  RW(  34  i , 
£IER,X,Y»M,PHI,3B(2»3),D(2,3J 

t.Xl,  X2,  Yl,Y2,Y3,Kl,K2.CC(2,3i,AAl2,3i  ,CON(5,5J 
INTEGER  C,1W(  13)  ,M1,M2, 1A,N,R 

£*****444*44*4 4  * *4* 4* ***** *4*4 4 *4 *4***44*444 4* 4  ********* 

C  VARIABLE  DEFINITIONS 

£*****444*  44**4  4 *4* 4** ** **  * ** ***** **** ******** 4 4* ** *** ** 

C 

AA.BB ,CC,0,U1 ,U2, VI , V2=AT  TRT  I  ON  COEFFICIENTS 
PSOL, QSCL  =  PRI  MARY  ANO  DUAL  SOLUTION 
M 1=NUME  ER  OF  ECUALITY  CONSTRAINTS 
M2=NUMEER  OF  INEQUALITY  CONSTRAINTS 
RW,IW=PARAMETERS  USED  BY  IMSL  ROUTINE  2X4 LP 
I  EP=ERROR  CODE  FROM  2X4LP 
N ® NU M E  E R  OF  UNKNOWNS 
A*MATRIX  DEFINED  IN  EQUATION  A. 3 
"LANCHESTER  EQUATIONS  ANO  GAME 
BY  F.H.  MOOSE  AN  C  J.  M  .  WOZ  ENCR  APT 
C  s  vECT  CR  CONTAINING  COEEFICIENIS 
OF  OBJECTIVE  FUNCTION  IN  LP 
Bs VECT OR  CONTAINING  THE  RIG 
OF  THE  CONSTRAINTS. 

I  A=M  1  +  M  2  +  2  ;I  S  THE  ROW  DIMENSION  Of  A  MATRIX 
QsCQLUMN  DIMENSION  OF  A 
Ml  *2 
M2=0 
N=3 

IA=M 1+M2  +  2 
Q  =  N+Mi+2 
P=M1 +M2 


C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 


OF  PAPER 
THEORY" 

OF 


[GHTHAND  SIDES 


c 

c 

c 


A  A  (1 
AAU 
AA  ( 1 
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non 


AA  (2 

,1) 

.60 

A  A  (2 

,2) 

= 

.900 

AA  (2 

,3) 

- 

0.0 

eea 

,1) 

s 

0.15 

66(1 

,2) 

9 

.10 

EB  ( 1 

,3) 

■= 

0.0 

EB  (2 

,1) 

s 

0.15 

65(2 

,2) 

S 

.3 

EB  (2 

,31 

* 

C.O 

CC(1 

,11 

s 

1.2 

CC  11 

,21 

s 

1.0 

CC(1 

,3) 

s 

15. 

CC  ( 2 

,1) 

s 

1.1 

CC  (2 

,21 

* 

0 .6 

CC  (2 

,21 

- 

15.0 

0(1, 

1) 

a* 

.00 

0(  1, 

21 

= 

.00 

0(1, 

31 

= 

0.00 

0(2, 

11 

s 

.00 

0(2, 

2) 

= 

.00 

0(2, 

31 

0.0 

C*****  **** **************** ********************* ********* 
C  REPETITION  LOOP  STARTS  HERE 

£*********************  ********************************** 
C 

7000  CCNTINUE 

ENTER  >  AND  Y 
WRITE  (6,29*1 

296  FORMAT  (IX, 'ENTER  X  AND  7  ONE  AT  A  TIME*] 

R  EAD  (  5  ,*  J  X,Y 


1000 

c 


61=  BElli  1  J*BB(2,  11 
62=  B  E  (1 ,2  )*B  6(2,21 
63=  B6  ( 1 , 3  1  +  8  B  (  2,  3  J 
01=  0(  1,1)*0(  1,21  *0(  1.3) 
02=  0  ( 2  » 1 1  +D  ( 2,21 *0(2,3) 


A ( 1, 1 Is  CC(1,  11 
/ X+1 5  « 

A ( 1 , 2 1 s  CC(1, 21 
/ X*1 5  . 

A  (1.31=  CC(1 ,31 
/ X+l 5  . 

A (2, 1 ) *  CC (2 , 11 
/ X+15. 

A  (2,21=  CC(2, 21 
/X+15  . 

A ( 2, 3  J  =  CC (2,  3) 


AAU,11  +  (Dl-Ull/Y  -  (Bl-Vll 
AA  ( 1 ,2 1  ♦  (Cl-Ull/Y  -  (B2-V2) 
AA  ( 1 ,3 1  *  (Dl-Ull/Y  -  (B2-V31 
AA ( 2 , 1 1  ♦  (C2-U21/Y  -  ( 8 1- V 1 J 
AA(2,21  +  (D2-U21/Y  -  (B2-V21 
AA  (2 ,3 1  ♦  (D2-U21/Y  -  (82-V31 


WRITE  (14.9591  ((I, J,A( I  1 ,J=1  ,2)  .1  =  1, 21 
FORMAT  (•  O'.  ,3X,'A(  •  ,13,13,  *1='  ,F  10 -31 

B  (11  =  1.0 
8(21*1.0 
CCNTINUE 

C  (11  =  1.0 
C  (21  =  1.0 
C (31=1.0 


130 


C 

200 

300 

C 


to  R  IT  E  (1^,1201 

FORMAT  (IX  ,  'CUEFFI  CIENTS  OF  CONSTRAINTS  (AIM 
DO  3  0  C  1*1,  R 

to  R  I  T  E  (  14,  2  00  1  (A <  I  ,K  1  ,K=1  ,N1 
T  h  I  J  NUMBER  =  N+Ml+2 
FORMAT! Ix, 7F15.U 
CONTINUE 

CALL  ZX4LP(A,  1A,B  »C,N,M1, M2,S ,PSuL ,DSOL,RW , 

I  W ,  i  E  F  J 


C  COMPUTE  Yl,Y2,Y3,Xl,X2,M 

C*****  *************  *******  **94  *******  *********  ********** 

c 

P  H 1  =  s 

Y  1=PSCL(11*Y/PHI 
Y2  =  PSCL  (2)*Y/ PHI 
Y3=PSCL(3> *Y / PHI 
X  l=QSCLi  1)*X/PHI 
X2=DSCl <2J*X/ PHI 
M  =  X*Y  / (PHI-151 
Kl=<  4  El-V  1J  *  Y1  * 

K2=(  ( C 1-Cli * XI  * 


( 82-V2 1 * Y2  ♦  (B3-V3J*Y3  *  M  1/X 
( U2-U2 1 *A2  -  M  i/Y 


Q************* 


COMPUTE  CONFLICT  MATRIX 


C* ******** **** 

c 


,»»*n***«<<  ****************************** 


2110 

2100 

C 


CONT 
CONT  INUE 


***************************************** 


AA(1,11*Y1  ♦  A  A ( 1 , 2  1  *Y  2  +  AA(1,31 


*  BB ( 1 , 1 i 
•»  BB  ( 1 , 2  1 


AA  ( 2 » 1 1  *Y  1  ♦  AA(2,2J*Y2  +  AA(2,3J 


i  11  *  Y1  ♦  0(2,11 

CC(1,  1J*X1  +  CC(2, 11*X2 


2)  *Y2  *  0(2,21 

CC  (  1, 21*X1  ♦  CC(2,21*X2 


*  0(2,31 


1,5 

,L  51  =  —  1  •  0*C(JN  ( L4 »  L5 1 


C  CN(  1 
*  Y3 

,11  = 

Cl  + 

COM  1 

,21  = 

0.0 

C  CN(  1 

♦  21  = 

AA  (  1 

C  ON  (  1 

,41  = 

AA  (  1 

CONI  1 

,51  = 

AA  (  1 

COM  2 

,11  = 

C.O 

COM2 
*  Y3 

,21  = 

U2  ♦ 

COM  2 

,21  = 

AA(  2 

com  2 

,41  = 

AA  (2 

C 0N(  2 

,5  1  = 

AA(  2 

C  ON  ( 3 

,11  = 

CC  (  1 

C  0N(  3 

,21  = 

CC(  2 

CON  (3 

,21  = 

VI  ♦ 

CON  (3 

,41  = 

C.O 

COM  3 

,5  1  = 

C.O 

CCN(4 

,11  = 

CC(  1 

C  ON  (4 

,21  = 

CC4  2 

COM  4 

,2  1  = 

C.O 

C  ON  (  4 

.4  1  = 

V2  + 

COM  4 

,51  = 

C.O 

C  ON  (  5 

*11  = 

CC(  1 

C  ON  (  5 

,21  = 

CC(  2 

COM  5 

,2  1  = 

0.0 

COM  5 

,41  = 

C.O 

COM  5. 51=  V3  ♦ 
DO  21  CC  L4  =  l,  5 

CO 

2  110 

L5  = 

CON ( L  4 
‘  INUE 


WRITE<14,1501 


li 


150 


160 


170 

C 

180 


190 

191 

19  2 


210 


C 

220 


C 

230 


240 


FORMAT (IX, 'RIGHT  HAND  SIDES  OF  CONSTRAINTS  (B)*) 
WRITE  <14, leOJ  ( B(  I  )  ,I  =  1,R ) 

_  ,  THIS  NUMBER  =M1+M2 

FORMAT  (1X»2F10.4) 

WRITE  (  14 ,170) 

WPITE(14, 180)  (C(I),I=1,N) 

FQRMAK1X, 'COEFFICIENTS  OF  OBJECTIVE  FONCT  ION  (C)  •  ) 

THIS  NUMBER  =N 

FORMAT  (1X.3F10.4) 

WRITEUP, ISO)  N 
WRITE  (1A, 151)  Ml 
WRITE (14 ,152)  M2 

FORMAT  ( IX  , 'NUMBER  CF  UNKNOWNS  (N)  =  •,120) 

FCRMAT  (IX,  'NUMEER  CF  INEQUALITY  CONSTRAINTS! Ml)  =  * 

>  I  5 ) 

FORMAT (IX, 'NUMBER  OF  EQUALITY  CCNSTRA I  NTS ( M2) =  • 
WRITE(14,210)  S 

FORMAT  (IX, 'VALUE  CF  OBJECTIVE  FUNCTION(S)  *' 

, F15 .6  ) 

WRITEUP, 2201  (  PSCL(I)  ,1=1, N) 

THIS  NUMBER  =  N 

FCRMAT  (IX, 'VALUE  CF  PRIMAL  SOLUTION  (PSOL)  =  • 

, 3F15.6) 

WRITE  (14,230)  (OSCLI  I )  ,  1=  1  ,R  ) 

THIS  NUMBER  =  Ml* M2 
FCRMAT  (IX, 'VALUE  CF  DUAL  SOLUTION  (DSOL)  »  • 

, 2F1 5  .6 ) 

WRIT  £(14, 2  AO)  I  ER 
FORMAT  (1X,'IER  =  •  ,15) 

WRITE(  14,570)  X 
WRITE  (14,580)  Y 
WRITE  (14 ,530)  XI 
WRITE (  14  ,540)  X2 
WRITE  (14,500)  Y 1 
WRITE (  14,510)  Y2 
W  R IT  E (  1 4 ,  520)  Y3 
WRITE!  14 ,550)  M 
WRITE  (  14 ,660)  K1 


570 

FCRMAT  (IX 

,  'X 

32 

1 

, FI  5.6) 

580 

FORMAT  (IX 

,'Y 

2 

« 

,F1 5 .6) 

500 

FCRMAT  (IX 

,'Y1 

s 

1 

, F15.B) 

510 

F  CRMAT  (IX 

,  'Y2 

2 

• 

,  FI  5.6) 

520 

FCRMAT  (IX 

,  'Y3 

s 

• 

,F15.6) 

530 

FCRMAT  (IX 

,'X1 

2 

1 

,F15.6) 

540 

FORMAT  (IX 

,  *X2 

22 

1 

,  F 1 5 . 6) 

550 

FCRMAT  (IX 

,  »M 

22 

• 

, FI  5.6) 

660 

FORMAT  (IX 

,'K1 

23 

1 

, FI  5.6) 

670 

FCRMAT ( IX 

,  *K2 

= 

« 

,F1 5.6) 

Q***** *************************************  ************* 

C  PRINT  PARAMETERS 

Q* ****** ******  ****************************************** 

C 

■  ‘  “  U1 

U2 
VI 
V2 
V3 


WPITE  ( 14,705) 
WRITE  (14,710) 
WPITE  (14,  7  15) 
WPITE(  14,720) 
WRITE (  14  ,  725) 
WPITE  ( 14 ,760) 
WRITE  (  14,750) 
w  R I T  E (  14,7  70) 
WPITE  (  14,750) 
WFITE(I4,7€D) 
WRITE  (  14,750) 
WRITE  (14,750) 


( ( A  A  ( 1  ,  J)  ,  J=1 , 3  ) ,  I  =  1 ,2 ) 
((BB(I,J),J=1,3),I=1,2) 
( (C  C  ( I  ,  J  )  ,  J  =  1 , 3  )  ,1  =  1  ,2) 
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70  5 

W 

FITE  ( 

14, 

7  50) 

(  (D 

1  I, 

J) 

t  J  — 

1,3) 

,1=1 

,  2 

F 

ORMAT 

(IX 

,'U1 

=  • 

,F1 

5. 

6) 

710 

F 

OR  MAT 

(IX 

,  'U2 

S  * 

,F  1 

5. 

6) 

715 

F 

ORMAT 

(IX 

,  'VI 

,  F 1 

5. 

6) 

72  0 

F 

ORMAT 

(IX 

,  *V2 

=  * 

,F1 

5. 

6) 

725 

F 

ORMAT 

(IX 

,  *  V3 

s  • 

•  F15. 

6) 

750 

F 

ORMAT 

(IX 

•  3F 1 

5.5) 

760 

F 

ORMAT 

(  IX 

,  *VA 

LUE 

OF 

A 

MAT 

FIX 

(AAI 

•  ) 

770 

F 

ORMAT 

(lx 

,  »VA 

LUE 

OF 

B 

MAT 

RIX 

<bB) 

•  ) 

780 

F 

ORMAT 

(IX 

,  'VA 

LUE 

OF 

C 

MAT 

FIX 

(CO 

•  ) 

790 

F 

ORMAT 

(IX 

,  'VA 

LUE 

OF 

0 

MAT 

RIX 

(Cl  • 

) 

c 

q***** 

c 

£***** 

c 


673 

674 
C 

£****♦ 

c 

c 


********************************** 
PRINT  CONFLICT  MATRIX 

♦  * *#4 ********  ******************** **** ********** 


WRITE  ( 14 ,6731 
WRITE  (14,6741 


urn  c  uiioiti  (  (CCM  L4.L3  j  ,L3=ii3  /  ,  i .‘♦•=1,31 
FORMAT  (IX, 'VALUE  OF  THE  CONFLICT  MATRIX  (CON)*) 
FORMAT  (lX»5F14.c) 

****************  ********************************** 
PRINT  MEANING  OF  I ER 

************************************************** 


»L  5  =  1 


,14=1,5) 


W  R  IT  E  (  1 4  *  3  SO) 

WFITE  < 14 ,4C0 ) 

WRITE  (  14,4101 
WFITE(  14,420) 

WRITE  1 14  ,4201 
WRITE  (  14,440) 

*  PIT  E  (  14 ,450) 

WRITE  (14,460) 
FORMAT  (  IX,  •  '  ) 
FORMAT  (IX, *IER=130 
F0RMAT(1X,'IER=131 

•  ) 

FORMAT  (IX,  *IER=132 
CONSTRAINTS'  1 
FORMAT  (IX, 'I  ER=  13  4 
UN6GLNC  EQ  * ) 

FORMAT  (IX,  *IER=135 
I NFEAS IELE'1 
F  ORMAT  (IX,  *IER=136 
SOLUT IONS') 

FORMAT  (9X,«00  NOT 


460 

C 

Q***** 

c 

Q***** 

C 


INCICATES 

INCICATES 


M2>N '  ) 
EXCESSIVE 


ITERATIONS 


INCICATES  REOUNCANC 1ES  IN 
INCICATES  OBJECTIVE  FUNCTION 
INDICATES  CONSTRAINTS 
INCICATES  PRIMARY  OR  DUAL 


SAT  ISFY 


CONSTRAINTS' 


**** ********************************************** 
REITERATION  ROUTINE 

********  *  ***************************************** 


9000 


WRITE  (6,470 
WRITE  (6,480) 

FORMAT  (IX, 'TO  CONTINUE  ENTER  l'l 
FORMAT  (IX, 'TO  STOP  ENTER  -I' 1 
R  EAD ( 5  »*  )  L7 

lF(ABS(L7*1.0).LT.1.0E-5)  GOTO  9000 
WPITE(  14,450) 

WRITE  (14,451) 

WRITE  (  14,452) 

FORMAT  ( IX  , '  •  ) 

F  ORMAT  ( IX , *******  *****  ***************  ******* i ) 
F  ORM A  7 ( 1 X , '  •) 

GOTO  700C 
CGNT  INUE 
STOP 
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APPENDIX  B 

BOBBIE  OF  EQUIL2BBIUH  SOIOTIONS  OF  THE  N*N  PEOB1EM 

We  first  consider  the  2*2  and  3*3  problem  and  extend  the 
results  to  the  N *N  problem. 

1  •  2*2  Problem 

In  general,  each  trivial  solution  leads  to  a  unigue 
equilibrium  solution  in  a  continuation  process  and  there 
will  be  equal  number  of  equilibrium  and  trivial  solutions. 
Section  4  discusses  what  happens  when  there  is  degeneracy. 
The  trivial  solutions  are  obtained  by  solving 


(eqn  B.l) 


Il(ValllS*a12Z4)  *  0 
22(u2*a2lVa2Zz4)  =  0 

W'hlV'zi^  ■  0 

24(u4*c12z1+c2222)  *  0 


At  first  glance,  there  would  seem  to  be  24  trivial 

solutions  corresponding  to  the  number  of  ways  of  maXing  the 

lefthand  sides  of  equation  B. 1  zero.  Each  lefthand  sides  can 

be  made  zero  by  either  maXing  z^  or  the  terms  in  parenthesis 

equal  to  zero.  But  closer  examination  reveals  only  six 

allowed  cases,  in  non-degenerate  cases,  corresponding  to  (1 
2 

*2  +  1)  =6  solutions.  Table  V  shows  how  these  cases 

arise. 
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TABLE  V 

Trivial  Solution  for  2*2  Probem 


Case  z.  which  are  made  zero 


Remarks 


1 

2 

3 

4 

5 

6 
7 


2 

2  =4  cases 


all  terms  in  parenthesis 
equal  to  0 
degenerate  case; 
see  section  4. 


2  •  3*3  prob  lem 

Here,  the  trivial  solutions  are  obtained  from 


T 

“1 

(U1 

+  a 

11: 

T 

'4 

+  a 

Z  2 

(u2 

+  a 

21  z4 

+  a 

Z3 

(u3 

+  a 

31 : 

4 

+  a 

Z4 

(u4 

+  c 

11: 

1 

+  c 

Z  5 

(U5 

+  c 

1 2 : 

1 

+  c 

*r 

“6 

(U6 

+  c 

13: 

f 

'1 

+  c 

12ValsV 

22z5*a25:6) 

32z5*a33:6) 

21z2*c31z3) 

22Z2+C52Z3^ 

23Z2+C35Z3) 


0 

0 

0 

0 

0 

0 


Again,  there  are  some  cases  which  are  not  allowed 
because  cf  inconsistencies  in  non-degenerate  cases.  Table  VI 
shows  the  different  cases. 

In  this  case,  the  total  number  of  allowed  cases  is 


l  1 


*  1  ] 


=  20 
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TABLE  VI 

Trivial  Solution  for  3::3  Problem 


Case 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 


z £  v;hich  are  made  zero  Remarks 


equals  to  0  • 

Z2’  Z3  ’  2 4  degenerate  case; 

see  section  4. 

Z1  degenerate  case; 

see  section  4. 
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j .  Extension  to  N*N  problem 

Novi  we  know  the  type  of  allowed  cases,  it  is  easy  to 
recognise  the  pattern  of  results.  The  pattern  looks  like  a 
Pascal  triargle  with  each  element  squared. 


N*N 

1*1 

l2  +  l2 

Number  of 
Equilibria  (N^) 

=  2 

2*2 

l2  +  22  +  l2 

=  6 

3*3 

l2  +  32  +  32  +  l2 

o 

II 

4*4 

l2  +  42  +  62  +  42  +  l2 

=  70 

5*5 

l2  +  52  +  102  +  102  +  S2  +  l2 

=  252 

For  N*N  in 

general. 

«.  ■i0‘ 

i-0 

N  can  also  be  written  as 
k 

■  ra 

The  proof 

of  the  identity. 

£(:)•■  6) 
i=n 

can  be  found  in  [Ref.  12] 


h .  Degeneracies 

Tor  the  2*2  problem ,  a  degeneracy  can  acise  when 
(or  any  ether  z  i)  is  made  equal  to  zero.  The  other  three 
lefthand  sides  in  equation  B.  1  are  made  egual  to  zero  making 
the  terns  in  the  parentheses  egual  to  zero. 

In  general,  such  a  case  does  not  correspond  to  a  new 
trivial  solution  because  there  is  an  inconsistency  when  z2  = 
-u3/c21  and  z2  =  “u4/c22*  However,  it  seems  that  by  making 
u 3/c2 1  =  u^/c 22*  the  inconsistency  no  longer  exists  and 
there  will  be  an  infirite  number  of  trivial  solution  as  long 
as  z3 ,  z4  are  chosen  to  satisfy  u2  ♦  a21z3  ♦  a22z4  *  0.  It 
turns  out  that  the  number  of  equilibrium  solutions  still 
remains  at  a  maximum  of  six  (disregarding  the  case  with 
infinite  number  of  equilibria).  Similarly  for  the  3*3  case, 
even  if  there  is  degeneracy,  the  number  of  equilibria  will 
not  exceed  20. 

that  the  above  is  true  can  be  shown  by  considering 
the  1*1  problem.  In  this  case,  the  two  non-degenerate  cases 
correspond  to  (a)  x  =  y  =0  (b)  v  ♦  cx  *  0  and  u  +  ay  =  0.  A 
degeneracy  can  occur  if  v  =  0  or  u  =  0  in  which  case  there 
seems  to  be  an  infinite  number  of  trivial  solutions  lying  on 
the  y  or  x  axis  respectively  and  hence  an  infinite  number  of 
equilibria.  But  when  the  actual  hyperbolas  are  plotted, 
there  are  only  two  intersections  and  hence  two  equilibrium 
points,  furthermore,  when  the  Continuation  method  is  used  to 
find  the  equilibria,  there  are  only  two  equilibria  irrespec¬ 
tive  of  whether  the  trivial  solutions  are  chosen  to  be 
degenerate  or  not. 


APPENDIX  C 

PROGRAM  LISTING  FUR  SOLVING  2*2 
USING  CONTINUATION  METHOD 


C 

c 

c  _ 

C  PROGRAM  LISTI  NG.  FUR. SOLVING_2*2  SYSTEM 

c 
c 
c 
c 

£******************  ************************************** 
C  THIS  PROGRAM  SOLVES  A  2*2  SYSTEM  FOR  THEIR 

C  ECUI LIBRIUM  POINTS  JSING  CONTINUATION  ME  T HODS .SEE 

C  SECTION  1  1 1 B  FUR  THEORY  AND  IMPLEMENTATIONS. 

C*****  4*******  *  ************************************  ****** 

c 

REAL  *  8  A.X.CEL  , T ,W K AR E A  .FNQRM , *K , U , R ,B  ,C»D#EP,AP, 
£>T,DX,OMAX,Xf  TEMP.TL  EFT  ,  T  $ ,F  1  ,F  2  ,F3  ,F4  ,E  RR 
PEAL  *  4  XA,XE,XC,XO,TA 

INTEGER  I  ,J,K,N,  M, MM, COUNT, I  A,I  DGT  ,  I  ER,NS  IG,  I  E,  I  REP 
£, CON, KK,NRECCN, SIGN 

£  ,FA$T,QCOUNT  , RUN , ERF  L AG , REFL AG, C HFL AG  ,CT R  ,CP1  ,JJ 
DIMENSION  WK  ARE  A  l  100  i  » WK. (  11 5  J  ,BP  (41  ,APt4,41, 
£XA(10G5)  ,  XE ( 1C05 J »XC i 10  05 } 

£,XD(  1C05J  ,TA<  100  51  ,DX<4J 
External  kn 

COMMON  /PEPLEN/  R(41 

COMMON  /PARAM/  A  (  2  ,2  J  ,B  (  2  ,2  i  ,C(2,2J,D(2,2),U(4) 
COMMON  /TREPAR/  XT  (4)  ,A(4J ,T  (41, 7  S,  FAST,  CON  ,  CGUNT  ,N 
£  ,OMAX,QCOUNT  ,  RUN  ,XFT  EMP  ,  SUM,  SIGN  ,  IER 
DATA  DEL/1. OD-16/, 

£MM/1/,M/1/,NS  IG/  5/  ,1  TMAX/200/,  I0GI/0/ ,  IA/4/  , 

£ERFL AG/O/  , REFLAG/O/ 

C 

C*** ******  *****  ****************************************** 

C  VARIABLE  CEF  INAT  IONS 

C *********  *****  *******  *********************************** 

C 

C  AP=S  EE  EQUATION  3-7  FCR  THE  MEANING. 

C  EP=AS  ABOVE 

C  >*THE  ROOT  OF  THE  2*2  SYSTEM  TO  BE  EVALUATED  EY  THE 

C  CONTINUATION  PROCESS. 

C  T*PARAMET  ER  DEFINING  THE  HUMOTOPYjHAS  VALUE  BETWEEN 

C  C  AND  I.C 

C  CEL=  SMALL  TIME  INTERVAL  USED  TO  APPROXIMATE  TFE 

C  PARTIAL  TIME  DERIVATIVE. 

C  AS,BS,CS  ,OS=THE  ORIGINAL  A,B,C,D  PARAMETERS  IN  THE 

C  LAfJCHESTER  ECUAT  l 'JN 

C  U,R=CCRRESFGND  TO  THE  SELF  ATTRITION  £  REFLENISH- 

C  MENT  CCEEF.  IN  THE  LANCHESTER  EQUATION. 

C  TS=P  ART  I T  ]0N  INTERVAL  FROM  T=0  TO  T=  1  IN  THE 

C  CONTINUATION  PROCESS. 

C  COUN T=COUNT£R  FOR  THE  *  OF  TIME  STEPS  ADVANCEC  IN 

C  THE  CONTINUATION  PROCESS. 

C  CON=  THE  CONDITION  FOR  REACHING  T  =  l 

C  FCN=  A  SUBROUTINE  USEC  BY  THE  IMSL  ROUTINE  ZSCNT. 

WK.NSIG  ,  ITMAX,IER,FNGRM=PARAMETERS  IN  THE  ROUTINE 

ZSCNT. 

IE,IA=PARAMETERS  IN  TFE  ROUTINE  LEQT1F. 

CTR=  COUN  T  ER  FOR  PLOTTING  THE  CURVES 
NR  EC  ON*  R  ERE  IITION  COUNTER. 


C 

C 

C 

c 

c 

C 


£***«* ********************  ******************************* 
C 

C  INPUT  ATTRITION  CCFF ICIENTS,  COUNTERS  AND  FLAGS. 

C 

Q********* *********************************************** 

A  ( 1,  1 )* 1 .000 
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c 

£*♦*** 
c 
c 

c***** 

c 


A(  If  2 1=  0 • 300 
A(2,  1  1  =  0 .6D0 
A(  2,  21  =  0. 900 
£11,  11  =  0.1 500 
£11. 21  =  0.  100 
£12,  11  =  0.  1500 
£12,  21=0 .300 
Ctlf  11  =  1  .2000 
Cll.  21=1.000 
Cl  2,  11=1.1000 
Cl  2, 21  =  0.600 
01 1,  11=0  .0000 
0(1, 21=0.000 
012, 11=0.0000 
0(2,  21  =  0 .000 
Ull)  =0.300 
U121=C.  300 
0(31=0. 100 
U(4)  =C.  200 
N=4 

T( 11 =0. 

11  21  =C. 

1(31=0. 

1141 =0. 

CO  UN 1=1 
CC0UNT=0 
RON=  1 
IER=  0 
IREP  =  1 
ARECQN=0 

********44*4*4****** 

CALCULATE  R  VECTOR  U 
SOLUTION  OR  MOOIFY  R 
*  * **  ***  *  *  *  *  **  ***  **** 


0  00 
000 
000 
0  00 


9998 

6 

C 

Q***** 

c 

c 

c 

c 

c 

£***** 

c 

7 

c 

£***** 

c 

£***«* 

C 


X  ( 11  =0. 6 153  6D0 
X  21  =0.  3  846D0 
X  (  31  =3. 0  "i  69DC 
X  (41  =0.9231  CO 
R(11=X(11*1U(11  +  A(1, 
C*X(3 1+B( 1  ,21 *X( 4 1 
F 1 21  *X(  2 1  *  ( U  (21  +  A  ( 2, 
£*X  (3  1+6(2, 21*  XI  41 
R(31=X(31*(U(i!+C(l, 
£*X(1  1+012,11*  XI 21 
R(41=X141*(U(41+C(1, 
£*X(1  1+0(2  ,21*X(21 
WRIT  Ell .99981 II, Xlll 

FORM  AT (  *  “  *  . . 

F12.41 


4  4*  ******$**$*>4444*44 

SING  ONE  POSITIVE  EQU 
VECTOR  AS  APPROPRIAT 

*****************4*** 


11*X(31+A(1,21*X(411+ 
11*X(3)+A12,21*X(411+ 
11*X(11+C(2, 11*X( 211+ 
21*X(11  +  C(2»21*X(211  + 
I.R(IJ . 1=1,41 


********** 

IL1BRIUM 

E. 

********** 


6(1,11 
6(2,  11 
Cl  1,  11 
0(1,21 


. •  1  •  K  l  1«  *  1  =  1  •  H  4 

0*,2X,*X(*,l3»*l  =  ',rl2.4,4X,  *R(*  ,13,'  1=' 


****************  *********************4************* 


WHILE  I  REF  (TFc  CCCE 
INUATICN  FRCCESS  WIL 
WILL  INCREASE  TS  IN 


FOR  REPEATINT  1  IS  1 
L  BE  REPEATED. EACH  RE 
MULTIPLES  OF  0.005. 


•THE  CONT- 
PE  TITIGN 


******* ****** 4* ******4* ******************  ********** 

CO  TO  8 


IF  (  •  NOT  •  (IREP.EQ.lll 
NRECCN=NFEC0N+1 

******************** 
SET  TRIVIAL  SOLUTION 

****  *********4****** 


************  ******************* 
************  ******************* 
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onooo  ooon  wonnoon 


X(1)=-C.10C0000/1 .20000 
X  (2J=-C. 000000/1.  10000 
X  <31  =-0.3000/  1.00000 
X  (4) =-C.CCC0/0. 30000 

t  ( n  =c  .oco 

T ( 21  =  C  .COO 
T  <31  =  0.000 
T ( 4) =  C . COO 
CHFLAG  =  C 
Q  C  0U  N  T  =  0 
R  0N=  1 
ERFL AG=C 
R  EFL  A  G  *0 
C  0  UN  T  =  1 
I  ER=0 
C  7R=  1 

TS=DFLCAT(NRECGN1 *0.005000 
C  QN=  IOiM  (  <1. OCO- T( 1 1  1/7S 1 +2 

************  4  4****  **********  **********************  ****** 

WHILE  T <  1  OR  CORRECTOR  STEP  DOES  NOT  RETURN  WITH 
ERRUR  FLAG  1ER 

*************  ********  *******  ********* *** ****** **** ****** 

IF(.NOT.((IER.NE. 129 1 .AND. < IE R .NE .130 1 .ANO . 

6  (  IER  .NE.131J  .AND.  (  COUNT.  L  E  .CON )  1 1  GO  TO  6 


TREAT  THE 
IS  SET. 


SINGULAR  CASE  IF  REFLAG 


ERFLAG 


CALL  TREAT (REF  LAG » ERFL AG»CHFLAG 1 


************************* 
PREDICTOR  STEP: 

SET  UP  THE  SET  OF  01 
CALL  LEGT1F  TU  SOLVE 
£*****  ******************** 
c 

DO  101  1=1, MM 
EP<11=<R<11 
EP  <  2  1=  (  R  (  2) 
eF<31=<  R<31 
EF  <4  )=(  R ( 4} 
AP<1  ,1 1  =U  ( 1 
AF< 1,21 =0.0 
AP<1  ,31  =  A(1 
AP  ( 1  ,4 1  =A  ( 1 
AP(2,1) =0.0 
AP(2,21=U(2 
AP<  2,3  1=A<2 
AP  <  2 ,41 =  A<  2 
AP  (3,11  =  C  ( 1 
AP(3,21=C(2 
AP ( 3  ,3 1  =U  (3 
AP l 3 ,41 =0.0 
AP ( 4 ,1 1 =C ( 1 
AP  ( 4  ,2  1  =C  (2 
Aft  4,31  =0.0 
AP (4 ,4) =U (4 
CALL  LECT1F 
CO  102  K=  1 , 
X  (K1  =  8P< 

102  CONTINUE 

101  CONTINUE 


******************************* 

FFERENTIAL  ECUATICN  AND 
THE  SYSTEM  AX=B 

******************************* 


102  CONTINUE 

101  CONTINUE 

C 

Q***** ********  ************ 

C  CORRECTOR  STEF: 


-8(1, 
-E<2, 
-D(  1, 
-D(  1, 
1  +A  (  I 
CO 

,11*X 

,21*X 

CO 

1 +A  ( 2 
,1)*X 
*  21  *X 
»11*X 
,11*X 
1+C(1 
CO 

,  2 1  *  X 

Sc1** 

1  ♦C  C 1 
<AP,M 
N 

K  1 +X  ( 


11*X<3)-6U,21*X(411*DEL 
ll*X131-B(2»2i*X(4)l*DEL 
l)*X<il-C<2,ll*X<211*jEL 
21*XI11-0(2,2)*X(21 1*DEL 
tl >  *X(  3  J ♦  A (  1 1 2)  *  X(  4 ) 

<11*8(1 ,1)*T(11 
(11+8(1 ,21*7(11 

,11*X(31+A(2,21*X(41 
(21+8(2  ,11*7(21 
(21+8(2 ,21*7(21 
(3  )  +0(  1,11*7  (3) 

(31+0(2  » 11 *T (31 
,11*X(1 1+C( 2 , 11 *X( 2 1 

(41+0(1  ,21*7(41 
(41 +0(2 ,21*T(41 

,2  1  *X(  1  1+C(2,21*X(21 
,N,  IA,dP,IOoT,HKAREA,IEl 


C*4444  4444 

c 


C  COMP 

C 


104 
C 
C 
C 

C 

C 

£444444444 

C  CCMP 

C********* 

c 


& 

£ 

£ 

£ 


1031 

C 

£*>44444444 

C  OUTP 

c********* 

c 


4  4444  4  4  4  4  44  4  4  4  44  4  44  4  4  44  4  444  *  444444444  4444444444 

CALL  ZSCNT  I FC  N  »NS  J  G»  N  »  ITM AX » T  ,X*FNCRM»WK,  IERJ 
OMAA=25C.QDO*DmIn1!DX1  1 J  ,  D  X  (  2)  ,  CX!  3  } ,  D  X  1 4)  ) 

UTE  DIFFERENCE  BETWEEN  PREVIOUS  AND  CURRENT  X 

DC  104  K=1  .N 

CX !  K  J-D  ABS!  XT! K )-X! K)) 

CCNTINUE 

CALL  SUBROUTINE  TO  DETECT  FAST  CHANGING  COMPONENT 
CALL  DETECTIDX , REFLAG, ERFLAG) 


****  *  4  *****  4  44  44  44444  4444*44  44444  4444  4444 

UTE  ERROR  IF  T  IS  NEAR  1.0 

4  444  44444444  4444  4  44  44  444  444  4  4  444  4  444  4  444  4  444  444 

IF!. NGT.(DABS!T!1)-1. 000). LE.O. 000500)  )  GO  TO 
1021 

f  1=-X!  1 )*!U ! 1)+A!1, 1)*X  <3)  +  A (1,2)*X<4)) 
■*Rd)“Bll»l  )*X!3)-B(1,2)*X!4) 

F2=-X!2 )*  IU ! 2) *A (2. 1)*X !3)+Ai2,  2)*X( A)) 
■*R!2)-B!2,1  )*X(3)-B!2,2  )*X!  4  ) 

F3=-X! 3  )*!U!3)+C<1  ,  1)*X  (l)K(2i  1)*X!  2)  J 
«R(3)-D!1,1  )*X!1)-D!2,1  )*X!  2) 

F4=-X<4 )*!U (4) +  Cd,  2)*X  <1J  +  C!2,  21*X(  2J  J 
♦  Rl4)-D(li2»*XU)-DI2,21*X(2l 
ERR=CSORT!F 1*42+F2**2*F 3**2*F44*2) 

WRITEU  ,997  mil.ERR 
CONTINUE 

*44*  444*4 44  4  4  4  4*  * 44  *44*44444  4***4 44*  *  4444444444 

UT  V £ LUE  CF  X  VECTOR  WHEN  T  IS  NEAR  1.0 

4  44  4 4 4444 444  44 44  4  44 4444  44444 444 4  4444  4 444  4 444444 


IF(.NGI.(T  di.GE.  (0.99 DO))  JGQ  TC 
“  ~  El  1 , 955 


865 

C 

C* 44 44  4444 
C  CHAN 

C4  44  4  4  4444 

c 


WRITEil ,955 )T !  1 ) 

WRITE! 1 ,999 ) !K,X!K) ,K=1 
CONTINUE 
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.N) 


6 
C 

C444444444 

C  COMM 

C+4444  4444 

c 

c 

c  c 

c 
c 
c 
c 


44444  44  44444  444444444444444444444444444444 44 4 44 

GE  DOUBLE  TO  SINGLE  PRECISION  FOR  PLOTTING 

44444 44 44 44 444 44 444444444444444444 44 4 44444444 44 

TA1CTR)=SNGL!T  (  1)  ) 

XA  (CTR) =SNGL (X  l  1) ) 

XB  !CJR)=SNGL!X  (2)  ) 

XC  <CTR)=SNGL!X  (3)  ) 

XD(CTR)  =SNGL!X!4)  ) 

COUNT  =  CCUN  T  +  l 
CTR*CTR*1 
GO  TO  5 
ONTINUE 

444444444 44 4444444444 44 44444  444 44  444  4 4444  444444 

ENCE  PLOTTING  ROUTINE 

444 44  444444 4  4444  4*44  44444444  44444  4 44  4  4444  444444 


ALL  TEK61S 

ALL  CONPRS 

ALL  FAGE!14.0,10.5) 

ALL  NOBRDR 

ALL  BLOWUP! 0 .4) 

ALL  AFEA2DU2.0,  8.0) 
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C  ALL  XNAME  ('  T  IME 
CALL  Y  NAN  E  (* 

CALL  XNONUM 
CALL  1NONUM 


SECS$',lOO) 


Xli)  $*,  100) 


INPUT  HEADING  IF  REQUIREC 


CALL  HEADIN!  •  INPUT  HEADING  HE  RE**** ******* ***♦*' 

, ICO,  1.  ,4) 

C  CALL  HEADIN!  •  ***********AND  H  ERE  **********  <****$' 

C  ,1 00 , 1  .  ,  A ) 

C  CALL  HEADIN!  *  ********* AND  HER £***♦ *** $  » , 100,  1 . ,■*) 

C  CALL  HEADIN! '  *********AND  HER  E****  **********$  ' 

C  ,100,1. ,4) 

CALL  GRAFIO.O ,' SC  ALE* ,  1.00,-56.0, 'SCALE* ,35.0) 
CP1=CTF-1 

CALL  C  UK  V  £  !  T  A  ,  XA,  CPI  ,  - 1 ) 

CALL  CURVE  <TA  ,XB,  CPI  ,-l) 

CALL  CURVE1TA ,XC, CPI ,-l) 

CALL  CUSVE!TA,X0,CP1,1) 

CALL  EN0PL10) 

WRITE  U  ,99101 

9910  FORMAT  ('O' ,3X  , 'ENTER  1  TO  CON  TINUE  '  ,3  X  ,  •  OTHER  NO. 
€  TO  STOP') 

R  E  AD  (5  » *  )  1REP 
GO  TO  1 
8  CGNTINCE 

CALL  CONEfL 
STOP 

955  FORMAT! 'O', 3X,'T  =  *  ,024.12  ) 

99  7  FORM  AT!  ■  0',3X,'T  =  '  ,0  24.  14.3X ,  'ERRCR=  '  ,D24.6) 

999  FORM  AT!'  O'  ,3X  ,'X  i«  ,13,'  )= *,024.14/) 

END 

£*****  *************************************************** 

C 

C  SUBPROGRAM  FCN  REQUIRED  EY  THE  ROUTINE  ZSCNT 

C- 

£******************************************************** 
SUBROUTINE  FCN!  X  ,F,N  ,  T) 

INTEGER  N 1 1  ,  J  ,NN 

REAL  *  8  ),F,R,U,A,B  ,C,D,T 

CIMENSION  X!4),F !4), T!4) 

COMMON  /REFLEN/  R(4) 

CCMMON  /PA RAM/  A12,2),B!2,2)  ,C(2  ,2)  ,C(2,2),U(4) 
F!1)=-X11)*(U(1)*A!1  ,1)*XI3)+A!1,2)*X!4)  ) +1 R  1  1) 

&-B  (1  »1)*>(3)-E(1,2I*X(4)  )*T!  1) 
F!2)=-X!2)*iU!2)+A(2,l)<‘X(3)+A!2,2)*X!4))4!R!2i 
C-B(2,l)  *X!3)-E!2  ,2)*X14)  )  *T  !  2) 

F(3)*-X!  3  )*1U13)«-C1  1  ,  1)*X!  1)  +C(  2 ,  1)*X  12)  J+tR!  3) 

C-Dll  ,1)*X(1)-0(2,1)*X(2))*T!3) 

F !  4)  =-X(  A)*!U14)  +C  (1  ,2)*X11)+C(2,2)*X!2)  )  ■»  l  R  (  A) 

L-Oli  ,2)*‘X1  1 ) -C(  2 , 2  )*  X  1 2  )  )*Tl  4) 

RETURN 

END 

C 

Q* ***************** ************************************** 

c 

C  SUBROUTINE  TO  TREAT  SINGULAR  CASE . WHEN  EITFER 

C  THE  “REVERSE"  FLAG, REFUG  CR  ThE  ERRUR  FLAG, 

C  ERFLAC-  IS  SET, THE  PARTITION  INTERVAL,  TS  WILL  BE 

C  INCREASED. IN  AUDIT  ION ,  IF  REFLAG  IS  SET, THE  FAST 

C  CHANGING  COMPONENT  WHICH  HAS  BEEN  FCUNO  NOT  TO 

C  FLIP  CORRECTLY  WILL  HAVE  ITS  SIGN  CHANGED. 

£********* *********************************************** 
C 

SUEROUT INE  TR EAT i P EF L AG, ERFLAG.CHFLAG  ) 


TO  CONTINUE', 3X  , 'OTHER 


c 

c 

c 

10 


c 

1011 


c 

c 

c 

11 


c 

1012 


c 

c 

c 


;2 


C 

1013 

20 


RE4L  *  6  X  7,  X  ,  T ,  T  S»0MAX,XF7£MP,SUM 

INTEGER  K.N, COUNT  .  CON  *  FAS T , QC CUNT ,RGN , SIGN  ,IER 

COMMON  / 7FEPAR/  XT ( 4  1  ,X l 4 1 , T ( 41  , T S ,F A ST , C ON , COUN T , N 
£  ,OMAX,QCOUNT,  RON, XFTEMP,  SUM,  SIGN,  IER 
DO  10S  K»1,N 

If  (.NGT.I  (  (REFLAG.NE.il  .AND. (  E  Rf  LAG  .NE  .  Ill 
£.0R. (CHf LAG.EQ.l 11 1  GC  TC  10 
7  (  K  1=T (  K 1  +  7  S 
XT  (  K1  =  X  (  K1 

IF  NEEO  TC  SET  X ( FAS  7 1=-XT( f AST 1  £  ERFLAG  IS  SET 

IF (.NOT.( (REFLAG. EC. 1J .AND . ( ERF  LAG  .  EQ. 1  ii i 
6  GO  7C  11 

>  (FAST1  =  -XT( FAST J 
X  7 1  FASTI  =  X (F  A  ST  I 

If  I. NOT.  (CHF  LAG. EQ. Oil  GO  TO  1011 
T  (Ki  =  T(K  J+TS 
T (21  =T ( 1  1 
T  (3 1  =T  ( 2  I 
T (41  =T ( 3  1 
TS=2  .000*  7S 
CHFL  AG-1 
T (K) =7 (K 1+TS 

CUN=  IDINT(( 1.0D0-T( 1 1 1/TS1  +  3 
CflUN  T  =  1 
CCN7INUE 
P  EfLAC=0 
E  RF  L  A  6=  0 

GO  70  2C 

IF  ONLY  ERFLAG  I  S  SE  7 

IF (  .NOT. ( ERFLAG .EQ. 11 1  GO  TO  12 
If  (  .NOT.  (CFFLAG.EC.011  GOTO  1012 
I(K1=T(K  )*TS 
T  (21  -T  ( 1  I 
T  (3 1  =T(2  1 
T (41 =T (3  1 
TS=2.000*7S 
CFFL  AG=1 
T(K1=T(K  1+TS 

C0N=I01NT((  1.0D0-T(  111/TS1+3 
COJN  T  =  1 
CCNTINUE 
REf  LAG=0 
EFf LAG=0 

IF  ONLY  REFLAG  I  S  SE  7 

GO  70  20 

IF(  .N07. (REFLAG. EQ. 11 1  GO  70  20 
X(FAST1=-XT(FAST1 
IF(.NQ7.  (CHFLAG.EQ.011  GO  TO  1013 
T ( Kl =  T(  Kl+TS 
T ( 21 =T( 11 
T( 3 1 =  7<  21 
T ( 4 J  *T ( 3 1 
TS=2 .OUO*TS 
CHFL  AG-1 
T(K1=T(K1+TS 

CON=  ID  1  NT ((  1 .000 -7 ( 11  1/TS1  +  3 
COUN  7  =  1 
CONTINUE 
REFL AG=0 
EPFL  AG=0 

CONTINUE 


94 


109 


CONTINUE 

RETURN 

END 


C 

C4‘4‘4<4*44*4‘4*<«44  4*4*44  4t  4  444  *****************  444  4  4*$*****$$ 

C 
C 

c 
c 
c 
c 

c*  44  44  44  44  4444  4  4  4  44  44  4  4444  44444 44 44 44 4444444444 44 44444444 

c 

SUBROLTINE  DETECT  4  CX,  REFLAG,  ERFLAG ) 

REAL  *  8  OMAX ,XFT EMP , SUM, XT , X ,CX( AJ.TS.T 
INTEGER  K  ,  CCU  UNT,  F  AS  T  ,  RON  ,  SIG  N ,  REF  LAG  ,  IER, ERFLAG 
N,  COUNT,  CCN 


SUBROUTINE  TO  DETECT  FAST  CHANGING  COMPONENT  OF 
X .UPON  RETLRNINGTC  MAIN  PROGRAM, THE  FLAGS  ERFLAG 
AND  R  EFLAGHI L  L  BE  SET  IF  ANY  SUCH  COMPONENTIS 
DETECTED 


C 

C 

c 

c 


L 

& 


,  II  f  LU  UIM  I  I  CUN 

COMMON  /TREPAR/  X  T  44  )  ,X(  4  J  ,T  (  4  J  ,  T  S  ,  FA  ST ,  CON,  COUNT 
, N,OMA>,GCUUN T, RON, XFTEMP , SUM, SIGN, IER 


IF  DXIKJ 
FAST  . 


DO 


>  OMAX  FOR  10  T S ,TH ATCOMPON  ENT  IS  CHANGING 


1015 

1014 

105 

C 

C 

C 


105  K  =  1,N 

If  (.NOT.  (DX(  K  i.GE.OMAXi  J  GO  TO  1014 
IF4.NQT.4 QCOUNT.GE.IO i )  GO  TO  1015 
FA  ST=K 

XF  TEMP  =  X4  FASTI 
RO  N  =  2 
QC  CUNT  =  0 
CONTI NUE 
(*C  CUN  T=QC  OUNT  ♦  1 
CONTINUE 
CONTINUE 

IF  X(FASTJ=-X7(FASTJ  MhEN  FLIPPING, SET  REFLAG 


10  21 

30 

10  22 

31 


40 
1017 
10  lb 
C 

r 

C 

c 


I F4.NCT. (RCN. 
SUM-XT ( FAS 
If  4. NOT. (4 
&  (XT4FASTIJ 

IF ( .NOT 
£  (XT ( FAS 

SIGN 
I  F(  . 

CCNT 
I  ER  = 
ERF  L 
GC  TO  4 
1 F  (  .  NO  T 
£  (  XT ( f AS 

S  1GN 
1F(  . 
R 

cont 

GO  TO  4 
1F4.  NOT 
I  ER  = 
ERFL 
CONTINU 
^ONTINUE 


NE.l 1)  GC  TO  101c 
T  i  +X  4  FA  ST  ) 

CABS ( X ( F AST )  )  .LE  .  ( 0 .  3C0#D  ABS 
J  J  .OR  .4  I  ER.GT.O  )  i  j  GC  TO  1017 
.4  4  DABS  4  X( FASTI l  .LE. 4 C.3DC*DAES 
Till  i.ANC.4  IER.GT.Oi  ii  GO  TO  20 
=  IDI  NT 4  XT (FAST) / X4 FASTI J 
NOT. (SIGN. LT. OJJ GO  TO  1021 
REFL  AG=1 
INGE 
0 

AG  =  1 
0 

.  (CABS4X4FASTI)  . LE. 4 0.3D0*DABS 

T)J)  ilGU  TO  31 

=101 NT4XV4FASTJ/X4FAST)) 

NOT.  (  SIGN. LT. OJl  GO  TO  1022 
EFLA  G= 1 
INUE 
0 

.(  IER.GT.OI  l  GO  TO  40 
0 

AG=1 

E 


CONTINUE 

IF  NO  COMPONENT 
FRCM  ZSCNT 


IS  CHANGING  FAST, BUT  THERE  IS  ERROR 
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LRNCHESTER'S  EOURTIONS  U.  .  (U>  NAVAL  POSTGRADUATE  SCHOOL 
HONTEREV  CA  A  B  NING  SEP  84 

UNCLASSIFIED  F/G  12/1  NL 


MICROCOPY  RESOLUTION  TEST  CHART 

NATIONAL  BUREAU  of  STANDARDS  1963-A 


10  20 


I F ( .NOT  « (  (RON . EQ. IJ.AND.I IEft.GT.QJJJ  GOTO  1020 
i  £R=0 
£FFLAG=1 
CONTINUE 
RETURN 
END 
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APPENDIX  D 

OSBIV ATIONS  OF  THE  HELATIONS  BETBEEN  ,  £„  AND  STABILITY 

*  V 

1 •  Neutral  Stability 


In  this  case,  £  =  e  =0  and  equation  4.5  reduces 

x  y 


to 


r,1(x‘xel)  +  (xy'xelyel^  +  =  0 

It  follows  that  both  sets  of  hyperbolas  merge  into  one  and 
all  the  points  on  the  common  hyperbola  are  equilibrium 
points.  The  lefthand  side  of  the  second  condition  in  equa¬ 
tion  4.4  can  be  manipulated  as  follows 

(u+aye) (v+cxe)  -  (b+ax£) (d+cye) 

=  ac  [  (nj+ye)  (d2+xe)  -  (W(n2+ye)] 

=  ac  ^xe(n1-n2)  +  ye(b2‘lJl)  +  niu2  '  n2yl] 

=  0 


This  result  implies  that  the  constant  term  of  the  character¬ 
istics  polynomial  D{s),  is  zero  ;  hence  one  of  the  eigenva¬ 
lues,  s1  equals  to  zero.  Factoring  out  the  characteristic 
polynomial,  we  get  the  other  eigenvalue  as 


For  points  on  the  first  quadrant  hyperbola, 


xe  >  -Ul  and  ye 

>  -r\  ;  therefore  s2  <0  and  neutral  stability  exists. 

Conversely,  s2  >  0  on  the  third  quadrant  hyperbola  which  is 
therefore  unstable.  The  results  are  summarised  as  fcllcws  : 

(1)  When  £x  ~  ey  ~  0/  inf initely-many  equilibria 
exist  as  points  cn  the  two  hyperbolas  on  which  x  =  y  = 
C; 

(2)  The  first  quadrant  hyperbola  is  neutrally  stable; 

(3)  The  third  guadrant  hyperbola  is  unstable. 

2.  Intersections  in  First  and  Third  Quadrant 

The  proofs  for  the  following  results  are  given  in 
this  section  : 

(1)  when  both  equilibrium  points  are  on  the  first 
quadrant  hyperbolas,  one  is  stable  and  the  other 
unstable 

(2)  when  one  equilibrium  point  is  on  the  first  quad¬ 
rant  hyperbola  and  the  other  on  the  third,  both  can  be 
unstable  or  one  will  be  stable  and  the  other  unstable. 

The  straight  lines  given  by  equation  4.7  are  plotted 
on  the  ex,  £y  plane  as  shown  in  Figure  D.l.  It  also  shews 
the  corresponding  regions  on  the  ex,  ey  plane  as  xg2  and  y^2 
vary. 

let  (x  ,  yel)  be  the  first  equilibrium  point  on  the 
first  guadrant  hyperbola.  The  first  stability  criterion  in 
equation  4.4  is  automatically  satisfied  since 

(u+aye)  +  (v+cxe) 

■  a(n1+yQ)  +  c(p,+x„) 


0 


Ml 


f°r  *xei*  yei^#  eguation  D.  2  represents  the  boundary 
line  of  stability  and  this  line  has  a  slope  between  n^/tx 
♦  Uj^)  and  (y gl  ♦  r|1)/*J1-  It  is  shown  in  Figure  D.2. 


/  s$*'/ 

/ 


/ 

/.<? 


"l 

/.•slope=  T^v 


^nl+yel^ 

&  /  yi 

^  / 


Figure  D.2  Boundary  line  of  Stability  on  Plane. 

To  investigate  the  stability  of  (xe2/  ye2),  we 
substitute  them  into  equation  D.2  and  obtain  the  condition 
for  stability  as 


ey  <U  e  W 

7  >  A 

s  (Pi+X  ,) 


(eqn  D.3) 


Zguations  D.2  and  D.3  are  only  different  in  the 
directions  cf  inequality.  This  means  that  on  one  side  of  the 
line  ey=  £x(n1+  Yei)/'{M1*  xei^  »  one  equilibrium  point  is 
stable  and  the  other  is  unstable.  On  the  other  side  of  the 
line,  the  opposite  conditions  exists. 

Note  that  (xe2,  YS2>  maY  or  maY  not  satisfy  the 
first  condition  of  equation  4.4.  This  condition  is 

a(ni+ye2)  +  c(ui+cx+xe2^  >  0  Oqn  D,4) 


From  Figure  D.1,  eguation  D.  4 
following  deductions  can  be  made 


and  earlier  results  the 


(1) 


ex>  °* 


x  „  >  y^0-  >  -n,  ;  hence  equation  D.u  is  satisfied. 


In  the  first  quadrant  of  the  ev,  e„  plane, 

is  sa 

in  Figure  D.3  is  stable  for 

plane  but 


*e2  '  "1'  Je2  '  "l  * 

The  region  labelled 


y  )  but  unstable  for 


(xe2'  ye2}  bat  unstable  f«r  (*el,  yel); 

(2)  In  the  third  quadrant  of  the  e  , 

x  y 

between  the  two  lines  where  x  _  >  0  and  y  _  >  -n,  , 

e  2  *  e  2  1  r 

equation  D.4  is  again  satisfied.  The  region  labelled 
is  also  stable  for  (X  „ , 

*  e2 

) ; 

plane, 

eguation  D.4  is  not  satis- 
so  both  (xel,  yel)  and  (xe2,  ye2)  are  unstable. 

in  Figure  D. 3 ; 


<x 

<3) 


el'  "el' 


Vt 


In  the  second  quadrant  of  the  e  ,e 
a  x  '  y 

y  <  -n  and  e  <  0, 

e  2  1  x 


Xe2  < 


fied 

The  region  is  labelled 
(4)  In  the  fourth  quadrant  of  the  e 


egua- 


e  plane, 

-•  y 

tion  D.3  is  not  satisfied  ;  so  (xe2,  ye2)  is  unstable 
tut  (xel,  y ei)  is  stable.  This  region  is  labelled 
mill  i-n  Figure  D.3; 

(5)  In  the  region  la  belled  ,  eguation  D.3  is 

not  satisfied  ;  so  (x  2*  Y-?)  is  unstable  tut  {xo 


y  )  is  stable. 

"  o  1 


"•el  ' 


tummmmrd 


v  el’  'cr 
unstable; 

(xe2.  ye2) 

stable 

as  above 

both 

unstable 

(xel’  >W 
stable; 

(xe2' 

unstable 


as  above 


vs 


OvVi) 

(ul*xel> 


Figure  D.3  Regions  in  e  ,  e  Plane. 


By  carefully  noting  the  signs  of  (xe2,  ye2)  and  the 
various  regions  in  Figure  D.3,  all  the  previous  deductions 
can  be  combined  ;  we  conclude  that  : 

(1)  When  both  eguilibrium  points  are  on  the  first 
guadrant  hyperbola,  one  is  stable  and  the  ether 
unstable; 

(2)  When  one  eguilibrium  point  is  on  the  first  guad¬ 
rant  hyperbola  and  the  other  on  the  third,  both  can  be 
unstable  or  one  will  be  stable  and  the  other  unstable. 
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3.  Two  Equilibria  in  the  Third  Quadrant 

The  simplest  way  to  show  that  both  equilitrium 
points  in  the  third  quadrant  are  unstable  is  t o  refer  to 
Figure  D.4.  Clearly,  we  must  have  x  <  -  v/c  and  y  <  -  u/a 
for  both  equilibrium  points.  In  that  case,  the  first 
stability  criterion  in  equation  4.4  is  not  satisfied  since 

(u+ayg)  +  (v+cxg)  <  0 
Therefore,  both  equilibria  are  unstable. 


-y 


Figure  D.  4  Two  equilibria  in  the  Third  Quadrant. 

4 •  Repeated  Equilibria 

This  case  corresponds  to  operating  exactly  on  the  r 

=  £x(yel  +  ni)/{*el  +  015  tbe  Ex'  Cy  plane*  The  Procf  is 

obtained  by  substituting  x  =  x  .  ,  y  .  =  y  into  equation 

4.6  and  solving  for  cy. 


(eqn  D.5) 


xe2  =  xel 


(Xei^i)  -  Vi; 


Rewriting  eguaticn  D.5,  we  have 


exyel  =  eyxel  +  ylEy  *  nlE: 


ex1=ey1+n-,e  -  u,  e 

y  el  x7el  lx  1  y 


subtracting  one  from  the  other. 


2exyel  =  2e-X-1  +  2Pl£-  ‘  2 ^ i E ' 


y  el 


ey  “  ex 


(eqn  D.6) 


Thus  we  have  shown  that  the  case  of  repeated  equi¬ 
libria  corresponds  to  points  on  the  straight  line  indicated 
in  Figure  D.3. 

Equation  D.6  can  be  also  obtained  by  setting  the 
second  stability  criterion  equal  to  zero  (  see  equation  D.2 
).  lhat  is  equilavent  to  saying  one  of  the  eigenvalues,  Sj 
is  equal  to  zero.  7he  sign  of  s9  is  determined  by  consid¬ 
ering  the  first  stability  criterion.  Following  the  same 
argument  as  in  neutrally  stable  case,  we  can  prove  that 
repeated  equilibria  on  the  first  quadrant  hyperbolas  are 
neutrally  stable.  On  the  other  hand,  repeated  equilibria  on 
the  third  guadrant  hyperbolas  are  unstable. 


APPENDIX  b 

BETWEEN  ST  ABILITY  AND  e  »  e  PLANE 

x  y 


VERIFICATIO NS 


Sene  representative  points  on  the  e  , 
chosen  and  the  corresponding  stablities  of 
calculated.  The  points  chosen  are  marked  F, 
Q,  T  and  W  in  Figure  E.l. 


plane  are 
the  equilibria 

G,  H ,  M,  N ,  P, 


Figure  E.l  Experimental  Verifications. 


i.  Procedure 


Having  chosen  the  points  on  the  plane,  we 

have  to  work  backward  to  obtain  a,  b,  c,  d,  u  and  v. 
Suitable  {*el,  yel)  ate  then  chosen,  followed  by  a  calcula¬ 
tion  cf  r,  s,  xg9  and  Y  e2-  txe2'  ? e2)  can  fce  stained 
directly  from  equation  4.6.  From  all  these  parameters,  the 
eigenvalues  can  be  calculated  and  results  compared  with 
theory. 


2 .  Fes  lilts 

The  results  are  tabulated  in  Table  VII.  They  agree 
with  the  theoretical  results  given  in  Figure  E.l. 


r\ 


TABLE  VII 

Results  of  Experimental  Verification 


APPENDIX  F 

PROG RAM  TO  PLOT  TRAJECTORIES  IN  1*1  SYSTEMS 


I MSL  ROUTINE 


C********* ************************************* #***<*****♦ 

C  PROGRAM  I  IP L  FORTRAN 

C  THIS  PROGRAM  PLOTS  THE  TRAJECTORIES  OF  A  1*1  SYSTEM. 

C  IT  CALLS  AN  IMSL  ROU TIN E  , DVERK  FOR  INTEGRATION. 

C  EEFOFE  EXECUTION , GET  ERMINE  THE  SCALE  OF  THE 

C  PLOT  AND  NUMBER  OF  CURVES  REQUIRED. 

C  THE  VALUE  .  "STEP"  IS  DETERMINED  EY  THE  RANGE  OVER 

C  WHICH  HE  RANT  TO  SET  THE  INITIAL  POINTS. CHECK  A, 6, 

C  C.D.AND  U<2).FIX  EITHER  X(1J  UR  Xt2i  AND  VARY  THE 

C  OTHEFjTHE  ONE  FIXED  HAS  TO  BE  RESET  AFTER  EACH 

C  CURVE  HAS  BEEN  DRAWN. IT  IS  RESET  IN  THE  LAST 

C  ASSIGNMENT  STATEMENTS  THE  DO  LOOP. TO  EXECUTE* 

C  ENTER  "  1 1  PL  "  .  THE  EXEC  FILE,"I1PL  EXEC" 

C  MUST  EE  IN  THE  DISK.1T  MUST  BE  EXECUTED  ON  A 
C  TERMINAL  ATTACHED  TO  THE  TEKTRONIX  616. 

C**** 1* ********  ******************************************** 

c 

C 

Q* **************************** **************************** 

C  VARIABLE  CEFIMTIGNS 

C**************  4  4  *******  *$:***  *****************  ************ 

C  >=  VECTOR  OF  LENGTH  2  CONTAINING  ThE  UNSTABLE 

C  EQUILIBRIUM  POINT. 

C  A,B,C,D* U*R*ATTRITIQN  CCEFFI ENT S . 

C  CAPX  ,  CAP  Y  *  CAP  IX  * CAP1 Y-ARRAYS  USED  TO  STORE  X"S 

C  FOR  PLOTTING  PURPOSES. 

C  £EL=  INTEGRATION  STEP  SIZE. 

C  IND,CC,W,TOL,TEND=PARAMETRS  REQUIRED  EY  IMSL  FOUTINE 

C  O  V ER  K 

C  1ER=  ERROR  MESSAGE  NUMBER  FRUM  DVERK 

C  OIR=£ONSTANT  FACTCR  DETERMINING  THE  DIRECTION 

C  OF  PERTURBATIONS  FROM  THE  UNSTABLE  POINT 

C  FCN1  =  EXT  ERNAL  FUNCTION  REQUIRED  EY  DVERK 

C  NSTEP=NUMBER  CF  CURVES  TC  PLCTJALSO  •  OF  STARTING 

C  STEP  =  INT  E  F  V  AL  BETWEEN  DIFFERENT  INITIAL  POINTS. 

C  INTEGRATION  POINTS. 

C  K*  KM  1— COUNTERS  FUR  PLOTTING  ROUTINE 

q********* *************************************  *********** 

C 

INTEGER  N.INO.NW , I ER ,K, NPOINT,KM 1 *KK  ,  JJ, K  1,K2  ,KM2 
INTEGER  N,INC»NW  ,IER  ,K . NPOi NT , KM1 ,KK ,  J,NS7£P 
PEAL  *  4  X(2),CCi24J  ,  W l 2  ,9) , T ,TOL . TEND, DEL, R(2i,A,B, 
&C,0,UI2i ,CAPX(1000J, CAP Y( lOOOI , X  A (  10 00 ) , X E( 10C0 J , STE P 
EXTERNAL  PCN1 
COMMON  R  ,  A  ,  E  ,C, D , U 

DATA  NH/i/tH/2/ , T/0. 0/,T0L/0 .001 0/ , I NO/17 , 1 ER/O/ 

C 

C  ENTER  ATTRITION  COEFFICIENTS  HERE 

C 

A=C. 6000 
£*0.  2000 
C=1.0000 
0*5.  CC 
U(li*0.60 
U(2J  *1.00000 
C 

C  ENTER  ONE  EQUILIBRIUM  POINT  HERE  FOR  THE  CALCLLAT ION 

C  OF  R  VECTOR 


X(1J*5.0CCCO 


ooooo  oooooooooooo  on  noooooooooo 


X(21=4.CCCCG 

=xu  |4(uui*a*x(  2)i+e*x(2i 

P(  21  *X121ML(2l  +  C*X(  lll+C*X(  1J 


ENTER  R  DIRECTLY  HERE  IF  THEY  ARE  KNOWN; REMOVE 
THE  COMMENT  CHARACTERS. 


P(l)=****< 

P(21«***** 


***** ******** * *******  ****  ******************************** 
SET  PLOTTING  PARAMETERS 

*****4*** **** 4  4*  ********** ** ********* ************* ******* 


999 


CALL  TEJU1E 

CALL  COMFRS 

CALL  VRSTEC(0.0,01 

CALL  PAGE  (l4.},ll.51 

CALL  NQBFCR 

CALL  CROSS 

CALL  BLOWUP  (0.51 

CALL  AR£A2D(13.7,9.501 

CALL  XNAME  ( 'TCTAL  X  FGRCES',1001 

CALL  YNAN  E  ( 1  TOTAL  Y  FORCES'  ,  1001 

CALL  FRAME 


INSERT  HEADINGS  WITHIN  CLOTESiREMOVE  COMMENT 
CHARACTER 


CALL  HEAOJN  (' *******************************$  • 
£,100,1. ,41 

CALL  HEAD  IN ( ' *****************************$• 
£.100,1.  ,41 

CALL  HEAD1NI •*********************$• , 100 


£.100,1. ,41 

CALL  HEA01NC ' *****************************$• 


£  ,  1 00 , 1 .  ,  4  1 

CALL  GRAF  (-3. 000, 'SC  ALE',  10. 0,-6. 00,  'SCALE',  1C. 001 


INSERT  INITIAL  CONDITIONS  FOR  INTEGRATION 


STEP  =  C« 5CCC 
NSTE  F  =  1 0 
XllJ  *1.0CCC 
DO  100  J*  1  ,NSTEP 

X  (21 =-5.00C+FLOAT(Jl ♦STEP 
K-l 

0  EL*C .01C 
NFOINTs-OO 

I  F  (  •  N  OT.  ( (IER.LE. 01. AND. ( IND . GE.O 1 .ANC.( K. IE. 

.  N POI N  Till  GO  TO  6 

TEND*FLCAT (K1*0EL 

CALL  OVERK  (N  ,F  CN1  ,T  ,X,  TEND  ,TOL,IND,CC,Nln,W,IER  1 
XA  (K  1=X (11 
XB  (K1*X(21 
WR  IT  E  (  fc  ,99  5 1  IE  R 


FOPF AT ( 'O'  ,3X, •  IER*', 131 


CAFX(K1*X(  11 
CAF Y(K1=X(  21 
K=  K  *  1 
GO  TO  5 
CONTINUE 
KM1  =  K-1 

CALL  CLRVE(CAPX,CAPY,KMi,31 


6 


on 


RESET  1NTEFGR  AT  ION  PARAMETERS  AFTER  EACH  CLRVE. 

T  =0.  C 
I  ER=0 
I  NO*  1 

X(1)*1.0C0C 
100  CONTINUE 

CALL  ENDFLIOJ 
CALL  OONEPL 
STOP 
ENO 
C 

£******* ******  * ***********  « ** **************  4* ************* 

C  SUBROUTINE  FCN1  RECUIRED  BY  OVERK 

£****« ******** 4 <*** ******* ******** ******** ***** *** ******** 

c 

SUBROUTINE  FCN1I  N,T»  X  »XPRIJ 
INTEGER  N#1.J 

REAL  4  A  X(N)  .XPRI(Ni,T,R(2i.A,BtC,0tU(2J 
COMMON  RfAtBtC.OiU 

XPRI  U)*-X(1J  *(U(  1  J  +  A*XI2  1  i*-R  til  —  E*Xl  2  ) 

XFRI ( 2I*“X(2i *(U( 2 ) *C*X( 1 li+R ( 2  J -C*X l li 
RETURN 

ENO 


c 

c 

c 

c 

c 

c 

£«**** 

c 

c 

c 

c 

c 

c 

c 

c 

c 

£****« 

c 

c 

c ***** 

c 

£***** 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C***** 

c 


APPENDIX  G 

FFGGFAM  TO  PLOT  BOUNDARY  CURVE 


******************** 

PROGRAM 
THIS  PROGRAM  CETERMI 
DOMAINS  OF  ATTKACTIO 
INTEGRATION  1NEGATIV 
FROM  AN  UNSTAELE  ECU 
BEFORE  EXECUTION  CHE 
TG  EXECUTE  ENTER  BKI 
MUST  EE  IN  THE  DISK. 
TERMINAL  ATTACHED  TO 
******************** 


*********************  *********** 

EK1  FORTRAN 

N ES  ThE  EOUNOARY  SEPERATING  THE 
N  IN  A  1*1  PROBLEM. BACKWARD 
E  TIME]  IS  USED  STARTING 
ILIBRIUM  POINT. 

CK  A,B,C»D,U(2J .AND  X(2J . 

.THE  EXEC  FILE,"EK1  EXEC* 

IT  MUST  EE  EXECUTED  ON  A 
THE  TEKTRONIX  618. 
******************************** 


*********  *********************** 
VARlAfiLE  DEFINITIONS 
******************************** 
X=VE CTOR  CF  LENGTH  2  CONTAINING 
EQUILIBRIUM  POINT. 

A.B»C,D,U,R  =  A7TRITION  CCEEFIENTS 
CAPX  fCAPY.CAPIX. CAP1Y- ARRAYS  USE 
FUR  PLOTTING  PURPOSES. 

CEL=  INTEGFAT ION  STEP  SIZE. 

INC, CC, W  ,TCL,TEND=PARAMETRS  REQU 
CVERK 

IER=  ERROR  MESSAGE  NUMeER  FROM  DV 
DIR=  CONST  AM  FACTOR  DETERMINING 
OF  PERTURBATIONS  FROM  THE  UNSIAB 
FCN1  =  EX  T  EFN AL  FUNCTION  REQUIRED 
E=CONSTANTS  EQUAL  1.0;FCR  SETTIN 
CF  PERTURBATION  TOGETHER  WITH  Dl 
K » K1 «K2  =  COUNT  ERS  FCR  PLOTTING  RO 
****  **************************** 


******************** 

******************** 
THE  UNSTAELE 


D  TO  STORE  X"S 


IRED  BY  I MSL  ROUTINE 
ERK 

ThE  DIRECTION 
LE  POINT 
BY  D V ERK 
G  OIRECT  ION  OF 
R 

U7INE 

******************** 


C 

c***** 

c 

£***«* 

c 


9998 

C 

c***** 


INTEGER  N.IND.NW  .  I  ER  ,K,  NPOINT.KM  1,KK  ,JJ,K1,K2  ,KM2 
REAL  *  <,  X(2i  ,CC  (241  ,«(2,9i  ,  T  ,TOL ,  TEND  ,U  EL,  E  (  2  J  ,  R ,  A,  B 
;,CtD»U.CAFMlCOOi»CAPY(lOOOi  » 
iCAPIXUOOC  J  ,CAP1  YUOOCI  ,CIR 
EXTERNAL  FCN1 
COMMON  R(2J 

COMMON  /PARA/  A,8,CfC,U(2) 

OATA  NW  /  2  /  ,N/2/  .  T/0.  0/ »  TOL/O  .001 0/ .  I ND/ 1  /  ,  IER/O/ 
i  .E/+1.000  , *1.030/ 

A=0.  7 
E=0.  A 
C  =  1.  0 
0=0.  6 
U(  1 J  =0. 15 
U(2J  =0.2 

►  ******* ********************************************* 
CALCULATE  R.S  VECTOR  USING  DATA  FR CM  LP  PROGRAM 

t ********  4  *******  ************************************ 

X(  1)  *£.70 
X  ( 2J  =2.00 

0(1)  •xm*(Ulll*A*XI2ll*E«(2I 
R C  2i -X4  2  J«  ( U<  2I+C*X(  imo*X(  U 
WRITE  (6, 9998  HI,  X{  U  , I , R ( I ) , 1=1 , N  J 

FORMAT!  *  0*  ,3X  ,'X  (  •  ,  I  3 » •  i  =  • ,F 12. A  ,4X  ,  •  R  I  •  ,  J3 , •  )  =  •  , F12  .  A 
*****************  ************************  ************ 
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c 

Q***** 

c 


INSERT  INITIAL  CONDITIONS  FOR  IN 

i  ********  *  ***********  ************ 

K  1=1 
K2*l 

CO  100  J  J*  1 , 2 


FOR 
R  ESP 


JJ«1  £  2, BACKWARD  INTEGRA 
ECTIVELY 


TEGRATIUN 

******************** 


TE  ALONG  ERANCh  #1,2 


C 

£***** 

C 

C 

c 

c 

£***** 

c 

5 

£ 


IFl.NOT.I  JJ.EG.1JJG0  TO  50 
D I P=0 . 02 

X(1J  =  X<1J*0IR*E(1  ) 

XI  2J  =  X(2l  +  DIR*El2i 
GO  TO  60 
CONTINUE 

0IR=-0.C3 

X(  1)=6.70«-CIR*E<1) 

X  1  2  1=2.  OOO+DIK  *EI  2  i 
CONTINUE 
K  =  1 
T  =  0.  0 
DEL=0 .100 

1FIJJ.EC.2)  0  EL=0  •  100 
NPOI NT=2G0 
I  N0=  1 

**************************************************** 


WHILE  NO  ERROR  ANC  NUMBER  OF 
LESS  THAN  NPOINT 


INTEGRATION  STEPS 


********  *  *******  ****************  *************  ******* 
I  FI. NOT. I  (IER  .LE.CI.AND.  I  IND.  GE.O  J  .AND.t  K.  LE. 


I  FI .NOT . ( IIER.LE.CJ.AND. I IND . GE .0 J .AN D.I K. LE. 

.  NPOINTJJJ  GO  TO  o 

TENO=-(FLOAT(K J*OEL* 

CALL  INTEGRATION  ROUTINE 

CALL  CV  ERK  (N,FCN1»T»X,  TEND  , TUL , IND , CC , N W, W , IER J 

PRINT  ERROR  MESSAGE  IN  ANY 

WR  IT  E  1 6  ,99  9  1  IE  P 

FORM AT  I *0'  ,3X, * IER  =  * , I  31 

STORE  X(1J  ,X(2I  FCR  THE  TWC  BRANCHES  OF  THE 
BOUNOARY  IN  TWC  ARRAYS 

IF  I.NOT.I  J  J.  EQ  .11  )GO  TC  501 
( APX1K) =X11 1 
(APYIKJ  =X(2  I 
K l  =  Kl+  1 
GO  TO  601 
CONTINLE 

CAP1XIK  )  =  X(  U 
(  AP 1 Y ( K  )  =XI  2J 
K2  =  K2»1 
CCNUNUE 
K*M1 
GO  TO  5 
CONTINUE 
CONTINUE 

SET  NUMBEF  OF  PLOTTED  PCINTS  TO  ONE  LESS  THAN 
THE  NUMBER  OF  DATA  POINTS 


ONE  LESS  THAN 


► 


c 


c 

Q4**** 

c 

£***** 


KM  1=  K 1-  1 
KM2=K2-1 

********  *4*******************************  ****4****** 

COMMENCE  PL CT TING  THE  BOUNDARY 

4*** ****4**** ******* *****************4**44********** 


8 


CALL  TEKtlfi 

CALL  COM F  F  S 

CALL  VRSTEZiO.O.  0) 

CALL  PAGE <14.7,1 1-5) 

CALL  NQBPCP 

CALL  CROSS 

CALL  BLOWUP  (0.551 

CALL  AREA2D1 13.7,9.0) 

CALL  XNAMEI* 

CALL  YNANE1* 


TOTAL  X  FORCES*. 100) 
TOTAL  Y  FORCES*  t 100 ) 


W 

8 

C 

c 

c 

c 

c 

c 

c 

c 


INSERT  HEADINGS  WITHIN  QUOT  ES  iREMCVE  COMMENT  CHARACTER 


C 

£***** 

C 

C*#*** 

c 


CALL 

»100 

CALL 

.100 

CALL 

.100 

CALL 

,100 

CALL 

CALL 

CALL 

CALL 

CALL 

STOP 

ENO 


FE A  DIN  1 ' *******************************$  • 
,1.,4) 

HEAD INI • *******  ********************  * 4$ • 

>1.  ,  ^  J 

HE A  DIN I  *  *******  *  ******* ******  100 

flat  4) 

HEAD  INI  *  *****************************$* 
f  la  «  A  J 

GRAF  16.000,  'SCALE* ,7.000,0.400, 'SCALE* ,4.000) 
CURViElCAPX,  CAPY  , KM  1,2) 

CURV  E1CAF1X.CAP 1Y,KM2,2) 

ENCFL(O) 

OONEFL 


********  ******************************************** 
SUBPROGRAM  CALL  BY  INTEGRATION  ROUT  I  N  E ,  D  V  ERK 
**************************************************** 

SUBROUT  IKE  FCN1  (  N,  T,  X  ,XPR  I) 

INTEGER  N ,  I ,  J 

REAL  *  4  X  IN)  ,XFRI(N),T,R,A,B,C,D,U 
COMMON  P  (  2  ) 

COMMON  /PARA/  A,B,C,D.U(2) 

XPRI  ( l)»-xll)  *1U1  1)*A*X12))*R(1)-E«x(2) 
XPRI12)=-X(2)*(U(  2JtC*Xll)i+R(2i-Dni  1) 

RETURN 

ENO 


.  i 

.  N 


>4 
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A' 


r- 


c 

c 

c 

c 

c 

c 

£44444 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

£44444 

c 

£44444 

c 

£44444 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C 

L 

c 

c 


APPEN01X  H 

PROGRAM  TO  OBTAIN  PAYOFF  MATRIX  AND 
OPT  I  MOM  M IX  ED  STRATEGIES 


4***  444444444  <4 4444444444444444444444444444444444444 

THIS  PROGRAM  COMPUTES  THE  PAYCFF  MATRIX  ANC 
OPTIMUM  MIXED  STRATEGIES  FOR  X  ANC  Y.  THE  ALGORITHM 
15  BASED  ON  THE  THEORY  GIVEN  IN  CHAPTER  5. 

THE  PAYOFF  FUNCTIONS  ARE  COMPUTED  FCR  60*60  X,Y 
BUT  ONLY  15*15  SAMPLED  MATRIX  IS  PRINTED. 

THE  CUTPUT  ALSO  INCLUCE5  MATRIX  FOR  LX.LY.AND 
FINTIM. 

BY  CHANGING  ALAMOA  TO  LESS  THAN  1 »ThE 
OEJECT1VE  FJNCTION  CAN  BE  MADE  TO  EMPHASIZE 
MORE  OF  FINTIM  IN  ACCORDANCE  WITH: 

A  I X,Y)  =  ALAMDA* (lY-LX)*(1-ALAMDAJ*FINTIM 
BEFORE  EXECUT  ICN,  CHECK  ATTRITION  CQEFFIENTS, 

QX.QY , /L AM CA ,  X IS, X2S, STEP  1, ST EP2, CARDS 
LABELLED  C  1,C  2  ANC  SCALE  IN  THE  PLOTTING 
ROUTINE. PERTURBAT IGNS  ARE  GIVEN  BY  CARCS 
LABELLED  C3  AND  CA. 

TO  EXECUTE, YOU  NEED  THE  EXEC  FILE  "LOSCP  EXEC" 

AND  A  TERMINAL  CONNECTED  WITH  A  TEK616. ENTER 
"LGSQf"  WHEN  REAOY. 

4  4  4  4  4  4  4  4  4  4  4  4  4  44  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  44  4  4  4  4  4  4  444 
4444  4  444  4  4  4  444444444  ******************************** 

VARIAELE  CEFINITICNS 

4444444444444444444444444444444444444444444444444444 

Xt2)=ARRAY  CONTAINING  X  AND  Y 
X1S.X2S=VARIBLES  USED  TO  SET  THE  FARTITONS 
OF  X  AND  Y  VALUES  TO  CALCULATE 
A(X,YJ. THERE  ARE  60  INTERVALS  FOR  X  AND 
Y  IN  THE  RANGE  0.2Q  AND  0.750. 

STEP1  ,  S7EP2-TO  OBTAIN  THE  CORRECT  STEP  INTERVAL 
AS  DESCRIBED  ABOVE 

XPRIN1  ,XPR  IN2  “VARIABLES  USED  TO  SAFEKEEP  INITIAL 
VALUES  OF  X  AND  Y 
TX,TY,P,S=AS  OEFINED  IN  CHAPTER  5 
RX,RY=(QX-*]  ,  4GY-Y1  RESPECTIVELY 
XL OSS  ,YLOSS=LX,LY  IN  CHAP.  5  „  , 

CC.NM ,W,TOL,T ENC, DEL,  INU= FARAMtTE P S  RECJIRED  BY 

IMS L  ROUTINE  ,DVERK. 

I  EC=ERROR  CODE  FR CM  OVEKK. 

Z MAT  =  MATRIX  USED  FCR  PLOTTING  SURFACE  A( a  «  Y i  . 

R ATLCS=FAYCFF  MATRIX 
F  INTIM=f  IN  ISH  TIME  MATRIX 
R  ATA  V  E= 1 5* 15  PAYOFF  MATRIX  FuR  OUTPUT 
XLGSPR=LX  MATRIX  FCR  CLTPUT 
YLOSPP’LY  MATRIX  FCR  CoTPUT 
FTPR=FINTIM  MATRIX  FCR  OUTPUT 
R OWM I N  =RL W  MINIMUM  IN  PAYCFF  MATRIX 
C  OLM  A  X  =  COL  UMN  MAXIMUM  IN  PAY  J  F  F  MATRIX 
AXMIN-MAXIMJM  OF  THE  THE  ROW  MINIMUM 
X  INMAX  =  MINIMUM  OF  THE  CULUMN  MAXIMUM 
RHS® ARRAY  CONTAINING  RHS  CF  CONSTRAINTS 
OECOEF*ARRAY  CONTAINING  COEhF ICIENTS  OF  OBJECTIVE 
FUNCTION  IN  LP. 

alamoa*nume-:r  BETWEEN  0  AND  1  TO  determine 
RELATIVE  EMPHASIS  CF  LOSS  AND  FINTIM 
IN  AIX.YJ 

SMMI  N  s  SMALLE  S  T  VALUE  OF  RCW  MINIMUM 
P  AY* 15*15  MATRIX  USEC  TO  SAFEKEEP  RATAVE  AND 
COMPUTE  THE  OPTIMIZED  STA  REG  I E  S 
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oooooooooooooooo 


AND  PLOTTING 
CEFINEU  AfcCVE 


NS  1,  NS  2SDI  MENS  ION  OF  PA YOFF  M  ATR IX  ,EU  UALS  AO. 

NCGO  =  FLAG  TO  INDICATE  END  OF  INTEGRATION. 

JJ,KK=  INDICES  FOR  MATRICES 
C  OUN  T  *  COUNTER  FCR  INTEGRATION 
K  KM1  =  NUMBE  R  OF  POINTS  PLOTTED 

F  CN1 >  2  1 3  ,4=SU 8RQO  T INES  RE  GDI  RED  BY  DVERKiTFEY 

CONTAIN  1*1  EJUATICNS  FCR  TFE  VARIOUS 
STAGES  FOR  THE  BATTLE 

T  *TI M  E 

XA,YA  ,ZA=ARRAYS  USED  FOR  STORING 
XFR1 N1 tX PRI N2  ,RA TLOS  AS 
ZCOUM*COUNT£  R  FOR  XA.YA 
JC=COUNTER  FOR  ZA 
T  IMS  1»T  1  IN  CHAP. 5 
T  1=T 2  IN  CHAP. 5 

NPLOT*NLMBER  GF  PLOTS  REULiRcO 

£*«***  4****+*#  4  4  *4*4#*  4***  *4*  *********  4***  ********  ******** 

REAL  X1S,X2S. X(2>  ,R,S,RX, RY,TX,TY  ,  U,  C  Y  ,  A  ,  C  ,U  ,  V , 

£  M (2*9  I  »T »TOL» TEND .DEL* XA( 36501 , YA < 365 C i , RATI  0 
£  ,  RATLCS(60,60  J , XP R INI t XPK IN2 , XT  1 , XI 2, R QTR AT 

£  *  ZMATUO  ,6Ci  .RORMIM  60  J ,  C  OLMA  X  ( 60  J  ,  AX F  IN,  X  INMAX 

REAL  F  IN  AL  «FI  NT  1M  ( 60 1 60i  .  T IME  1  «RAT  AVE  (  20 , 20  I 
£  , SMMI N,PAY(17 ,32J  , JEMPG,RHS( 1 5 J ,UECUEF (15 i  ,AL AMDA 

£  , XLOS (6C.6C) , YuCS (60 , 601 ,FTPR (20,201  , 

£  Y LOS P P  (20,20 J  ,CC( 24 J ,T1,STEP1 , ST EP2 , Y LCSS , XLOSS 
REAL  YLCXL,ZA(3650J ,XLCSPR(20 ,20 J 
INTEGER  INO, NW,  IE  D,K,CGUNT,NSTEP,NS1  j  NS2,  J.NOGO 
£  ,  ZCULIN1 ,  NX  ,NY  ,NPU  I  NT  ,  J  J  ,KK,  JN  ,  JC  ,  KKMl 

COMMON  /PARA M/  A,  C  ,U  ,  V  ,S  ,  R  ,B  ,  U 
EXTERNAL  FCN1 
EXTERNAL  FCN2 
EXTERNAL  FCN3 
EXTERNAL  FCN4 
C 

8************* * 4  **************  *  ********************  ******* 

SET  PARAMETERS  ANO  COUNTERS 

£#*****  +  ♦*****<  *4  *********************************  ******** 

C 

A=0. 70 
C*1.00 
U=0. 15 
V  =  0.  200 
8=0.4000 

o«C.t>ao 

A  L AMO A  =  1.00 
F  INAL  =  0.10 
N  W=2 
N  =  2 
T  =  0. 0 
TCL=0  .001 
I  NC=  1 
DEL=O.CI 
C  0  UN  T  =  1 
QX=10.CC 
Q Y  =  7  .  C  C 
X1S=3.C 
X  2  S=  3  .  C  0 
JC=0 
N  QGO=  C 
NSl=cC 
NS2=6C 
N  S=N  S  1  -  I 
Z  C  OU  N  T  =  Q 
NPOI  NT  =  NS1*NS2 
C 

Q***** *****************************************  *********** 

C  START  COMPUTING  &C*60  PAYOFF  FUNCTIONS 


£#**** ********  * 

c 

DC  100 
STE 
DO 


CARD  Cl: 
CARD  C 2 : 


CARO  C3 : 
CARO  C4: 


£************** 

C  STAGE 

C* ************* 

c 


|  4 4*4  4** *********************  4 *** *********** 

3  J*1,NS1 

:F1=(  FLOAT  (J-NSlli/39.33 
101  K=1,NS2 
JC=JC+1 
T«O.C 
COUNT-1 
IN  0=  1 

STEP2=( FLCAT(K-NS21 1/54.428 

X(11=X1S*(2.333-*STEP11 

X  ( 2 1=X2  S*(1.S17*STEP21 
XFRIN1  =  X(1J 
X  PR IN2=  X ( 21 

X  A  (  J  Cl  =  XPRI  N 1 
YA ( JC1  =  XPRI N2 

R  =  X(  U*(U+A*X(21)+B*X(21 
$  =  X(2)*(V+C*X(lli+D«X(ll 
RX*QX-X(  11 
R*  =  QY-X(  21 

XI  1J=XPRIN1-0.01 

XI  21-XPRIN2-0.05 
T  X=RX/R 
T  Y*RY/S 

1 *********** *** ********* ************* ******* 
1 

>4***4**  *********************4************** 


5000 


1028 


TENOO.  0 
TIME1=AMIN1  ITX 
1FI  .NOT.!  ITEND 
GO  TO  6000 

TEND*FLQ AT  I 
CALL  OVERKI 
.Nrif  Wt 1ED1 
IF(  .NOT.  I  IX 
.LE.  I  FINAL* 
1FI.N0T. 
XLGSS 

y  l  oss 

GO  TO  62 
CONTINUE 
XL  OSS 

yccss 

CONTINUE 
NLiGQ=  1 
XLCSI  J  ,K 
YLCSl  J.K 
RATLOSI  J 
F  INTI  Mi  J 
IF!  .NOT. 
SIGN* 
GO  TO 
CONTINU 
SI  GN= 
CCNT  INU 
RATLOSIJ 
-  A  L  AM  C  A 1 
CQNT INUE 
JgJN  T= COUNT 

GO  TO  5000 


#TY1 

. LE.T I ME  11. AND. (NO GO. EC. Oil) 
CCUNT  1*DEL 

NtFCN4tTiX,TENDiT0l,IND»CC 

1 1 1. L E. I F  INAL* CX )  1. OR.  (XI 21 
CY111 1  GO  TO  1028 
I  X(l) . LE . (FINAL *UX11 ) GC  TO  52 
*( 1.-FINAL1*UX 
*T*S+XPRIN2-X (21 


=7  *R*XPRI N 1-X (11 
*( l.-FINALl*UY 


1*XL0SS 

1*YL0SS 

,K1=YLCSS-XL0SS 

»K1*T 

( RATLGS ( J  .Kl.LT.O.QilGQ  TO  153 
1.0 
163 
E 

-1.0 

E 

»K )  =  A  LAMDA+RA7LCS I J.K 1 *11. 
*SIGN*FINTIM( J ,K1 


non 


CGNIINUE 
CALCULATE  T1 


6000 


50 

60 

***** 


w 


IF  (  .NOT  •  ( TX. LT  •  T Y J  J  GO  TO  50 
T  1=  T  Y- TX 
GO  TC  60 
CONTINUE 
T  1=T  X-TY 
CONTINUE 


C 

c«**** 

c 


******** ******************************************** 
STAGE  2 

****************  ************************************ 


L 

L 


154 

164 


1011 

51 


L 

L 


155 
16  5 


TENC*0.0 
7*0. C 
CQUNT=1 
INO*  1 

IF(  .N07.I (TEN0.LE.T1).AND.< 
GO  TO  6 

TENO*FLQAT(CCLNT  i*DEL 
If(  .NOT.  (  7X.LT.TY1  JGC  TO 
CALL  DVERK (N.  FCN1  ,T,X 
fCCtNMtW.IEOl 
IFUNOT.UXm.LE  .(FI 
(X(2J  .LE.(FINAL*QY)J) 
XL0SS  =  CX-X(1) 
YL0SS=(T*S)*XPRIN2 
RA  TLOSl J»  K I =YL  CSS- 
XL  0  S  (  J  ,KJ  =  XLOSS 
YLOS(J»KI=YLOSS 
FINTIMI J,K)=TX*7 
IF  (  .NOT  .  (  RATLO  S  (  J  > 
154 

SIGN*1.0 
GO  TO  164 
CONTINUE 
SIGN*- 1.0 
CONTINUE 

RA  TLOS ( J*  K  )= ALAMO A 
( 1 .-ALAMDA I*SI GN*F 
NO  GO* 1 
CONTINUE 
GO  TO  61 
CONTINUE 

CALL  DVERK4N.FCN2 »TtX 
.  CCfNM.W.  I  EL>  J 
IF(  .NOT.  I(X(1).LE.(FI 
1 XI2J .LE.(FINAL*OYIi i 
XL  OSS* ( T*R I+XPftiNl 
YL  OSS*  C Y-X ( 2J 
RATLOS(JtKi*YLOSS- 
XLOS(  J,KJ  =  XLOSS 
YLOSi  J,KI  =  YLOSS 
FINTIMiJ.K)*TY+T 
IF (.NCT.(RATLOS(J, 
TO  15  5 

Si GN*1 . 0 
GO  TO  165 
CONTINUE 
SI GN=-1  .0 
CONTINUE 

RATLOS ( Ji K i*AL  AMLA 
1 1 .-ALAMDA i*SI GN*F 
NO  GO* 1 
CONTI  NUE 


NOGU.EQ.OTJ  ) 


51 

t  TENC  fTOL  •  INO 


N AL*C X)  )  .OR. 
1  GU  TO  1C11 


->I21*(TX*SJ 

X  LOS  S 


KJ.LT.O.OJJGO  TO 


•  PATLCSI  J  »K J -*• 
INTIM  |J»KJ 


t  TENC fTOL  » I  NO 


NAL*  CXi  I  .OR. 
)  GO  TO  1013 
-xd  i*ar*Ri 


XLOSS 


KI.LT.0.0  JiGO 


♦  RATLOSI  J»K1* 
INTI M(J,K1 


1013 


6 

C 

£***** 

C 

0***4 

c 


500 

1005 


1014 


600 

C 

O**** 

c 

0***4 

c 


101 

100 

c 

£4***4 

c 

0**4* 

c 


1017 

302 


CONT INUE 
COUNT«CUUNT  *1 
I  ED=  0 
GO  TO  5 
CGNTINUE 

*********  4*4*4**  4*4* *** 4 ******** *4* 44 *4* 4 *44 4*4444* 

STAGE  2 

*444*4*4  *  *  ******  4*4*444*44***44*  44*4444*4  *44*444*4* 

1*0.0 
C0UNT=1 
1ND=  1 

If  (.NOT . (NOGO.EG.OJ JGQ  TO  6C0 
CONT INUE 

TEND=f  LOAT  (COUNT  i  *CEL 

CALL  CVERK(N,FCN3 ,lfX, TENC.TOL  ,IND 

•  CC  ,NW,W  ,  IEDJ 

1F(.NQ).((X41).LE.(FINAL*CX)J.CR. 

( X( 2) .Le.(F INAL*UYi)j )  GO  TO  1014 
XLGSS-CX-X(11 
YL  C  SS  =  C  Y-X 12) 

RATLOS( J,K)=YLCSS-XLUSS 

XL  QS( J  » KJ  =  XLUS S 

YLOStJiKj  -  YLOS  S 

FI NT1M(  J.K)=AMAXl ( JX,TY )  +  T 

IF  (  .NOT.(RATLOS(  J  ,  K  )  .L  T  .0 . 0  J  J  GO 

TO  156 

SI  GN=1  •  0 
GO  TO  166 
CONTINUE 
SIGN=-1.0 
CONTINUE 

RA TLUS ( J t  K )=AL  AMDA4RAT  LCS ( J  .K  J ♦ 

( 1 .-ALAMDA )*SI GN*F INT I M ( J  »  K J 
NOG  0=1 
CONTI  NUE 
C  GUNT  =  CG  UNT+ 1 
I  ED=0 

I  FI.  NOT.  ((X(li.LE.(FlNAL*CX)  j.OR. 

(X(2i. LE. (FINAL*QY>  UJ  GO  TO  1005 
GC  TO  500 
CCNTINUE 

****4**4  4  *44*44*  ****************44**444*4***444**** 

END  OF  EATTLE 

44**  44*444  44****  *****  4*  44  *******  4  44*4  *44**444  *4*44* 

XA(2CQUNn=XPRINl 
YA ( 2C  QUNT ] =XPR I N2 
NOGO'O 

Z  A  ( JC  i  =  RAT  LQS( J  *K  ) 

C0N1INUE 

CONTINUE 

44**44*4  4  *4 4* *******  44**4***4***4444444*4  4**4  4*444* 
COMPUTE  ROW  MINIMUMAND  SMALLEST  RCWMIN 

4444  4444  4  4  444  44*  4444444444*4444*4444444444*4*444**4 

DC  301  KKM.NS1 

RCMMIN(KK)  =R AT  LOS  ( KK»  1  i 
DO  302  J  J=2*  NS2 

If (.NOT .(RATLOS(Kk»JJJ .LT.RCWMINIKKJ i J 
GO  TO  1017 

RCWMINIKK j=RATLOS(KK ,JJJ 
CONTINUE 
CONTINUE 


CONTINUE 
DO  305 


lmj  jus  K  R-  2  » i\  a  i 

IF4.NGT  .(ROW MI  N(KK)  .LT.SMM 
SMMIN=RCWMI  N(NK) 

C  1NSER=KK 

1019  CONTINUE 

305  CONTINUE 

£************* 1***4  ****************** 

C  SAMPLE  15*15  MATRICES  FROM  601 

£***** ********  *  ******* ****************j 

C 

CO  301  JJ»l, 15 
J»A* J J 

00  3Cti  KK=1,15 
M4*Kk 

RAT  A  VE(  JJ„KK  )-RATLOS(J,  J 
PAY(JJ,KK)=RATAVE(JJ»KK, 
*LOSRK(  JJ,KK)  =  XLCSl  J»KJ 
YLOSPKI  JJ  ,KK)  =  YLCS(  J,K) 
f  TPRI JJ  ,  KK)  =F1  NT IM( J  ,K) 
308  CONTINUE 

307  CONTINUE 

C 

£***«* ********  *  * ***  *************** 4***1 

C  OUTPUT  MATRICES 

(;«****  *********  **************** *******! 

c 

WRIT  E (  1  ,580 (  (RATAVEI  JJ,KK)  ,M 
WRITE  (  1,58141 

WRITE  ( 1  ,58111  (  <XLOSPR(JJ,KK)  ,1 
WRITE(1,5614) 

WRITEI1, 98121  (  (YLOSPR  (JJ.KKJ  ,t 
WRITE(  1.9814J 

WRITE!  1,98  13 1  ((FTFR(JJ,KK),KK- 
WRITE (  1 ,9814) 

980  FORMAT!' 1' ,15F7.3//) 

9811  FORMAT  ('1'  ,15F7.3//i 

9812  FORMATCl' ,15F7.3//) 

9813  FORMAT  ('I*  .15F7.3//1 

9814  FORMAT  CO*  ,» - •) 

C 

£***** ********  ********  ****************1 

C  CALCULATE  OP  SOLUTION 

C  ACO  AES(SMMIN)  TO  ALL  -PAY(J,J 

C  POSITIVE 

£*** **********************************  i 

SMMIN  =  4BS(SMMIN)  +  2. 

DO  405  J  J*  1, 1  5 
RHS ( J J ) =1 . 

08CCEF(JJ)=1. 

00  «C4  KK* 1,15 

F AY ( JJ , KK 1= ( PAY ( JJ, KKl* : 
404  CONTINUE 

403  CONTINUE 

C 

C  CALL  SU8R0UTINE  LF  TO  COMPUTE 

C  STRATEGIES 

C 

CALL  LP(PAY,RHS,OBCOEF  ,SMM1N) 


IN) 1  GO  TO  1015 


******************* 
*60  MATH  1 4 
******************* 


*************  ****** 
******************* 
K=1,15),J«;=1,  15) 
KK=1 ,15)  ,-J  =  l  ,15) 
KK- 1,15) ,JJS1«15) 

- 1 , 1 5  )  ,  J  J  *1 , 1 5) 


******************* 


K )  TO  MAKE  GAME  VALUE 

******************* 


SMMIN) 


OPTIMIZED 


r********* ******  ****** *********************************** 

C  PLOTTING  ROUTINE 

C***«* **********  ********************** 

C 

CALL  TEK61E 


******************* 


ononoonon 


987 

986 

C307 


N  PLOT  *  1 

DO  102  L*1,NPL0T 

CALL  FACE (8. 50 ( 11*0) 

CALL  NCEKOR 
CALL  ELCWUPU.O) 

CALL  AREA2 Dt5.50,7.00) 

REACTS, *)XVU 
WRITEL6,906)XVU 
REACL5 ,*) Y VU 
WR  IT  E  (  6  »90  7 }  YV  U 
REAC  <5,*i ZVU 
WRIT£L6,908)ZVU 
CALL  FRAME 
CALL  S CMP LX 
CALL  XAXANGI45.CI 
CALL  YAXAN  G(  45 .0) 

CALL  ZAxAN G ( 45 • 0) 

CALL  X3NAMEC • X-CEPLOYMENT$ *  ,100) 

CALL  Y 3 NA M E ( 'Y-CEPLQYMENT* • , IOC ) 

CALL  Z3NAMEL  'PAYOFF  TO  X$«,100) 

INSERT  FEAOING  IN  4******  » REMOVE  "C" 


CALL 
CALL 
CALL 
CALL 
CALL 
7.01 
CALL 
CALL 
CALL 
5. 75 
CALL 
CALL 
CALL 
CALL 


HE ADI N( ************** $• 
hEADI N ( ************** $• 
hEADlNi **************$* 
FEADI N  ( '  *44*4*4** ****$• 


, ICO, 1 .C|4) 
,  100,1.0,4) 
,100,1.0,4) 
100,1.0,4) 
MESS  AG  ( •  *> 4 ***********  $  *  ,  100 , 1 . 2, 


102 

CONT INUE 

CALL  OCNEPL 

906 

FORM  A7(*0*,3X,'X  VU- '  ) 

907 

FORMAT  ('O',  3X,'Y  VL=*  ) 

908 

FORMAT  ( 'O'  ,3X,  'Z  VL=*  ) 

C 

WRITE  (  3  ,586)  1  NSER,  INSEC 

C988 

FORMAT  L*0',3X,'INSER=' 

VGLM3 Oil. ,1.,1.) 

VUAdS(XVU,YVL,ZVU) 

GRAF3  012.  50, 'SCALE',  7. 000, 2. 50,' SCALE*  , 
,-5.0,  'SCALE*  ,5.000) 

BCX3D 

RASPL N ( 0.  ) 

SURMA T(RATL0S,1,NS1, 1 ,NS2  ,0) 

ENCPL (0) 


WRITE  i  1 ,567)A  ,B,C  ,C 
FORMAT  CO'  ,2X  -• 


14,'  INSEC=* ,14) 


,-A,E,C»C  =  ',  4F13  .3) 


WRITE(l,S8fc)U,  V  ,a  L  AMOA 
FGPMAT('0',2X  ,'U,  V  =  '  ,2F13.3,2X,  • L AMDA* • , F 1 3.2) 
CONTINUE 
S  TCP 
END 


^***** ********  ******************************************* 
C  SUERCL  TINE  FCN1 

Q********* *********************************************** 

c 

SUEROUTINE  FCN1 (N  ,T,X ,XPR IME) 

INTEGER  N 

REAL  X  IN  )  »  XPR  IMEI  N),  T 
COMMON  /PAPAM/  A,  C  ,U  ,V  ,S  ,R,B  ,  0 
XPRI M  E  1 1 ) «-X ( l ) *(  U*A*X 12 ) )-B*X<2) 
XFRIM£I2J*-X(2)*(  V*C*XUI  )+S-C*X(l) 

RETURN 

END 

C 

Q***** ****************  *********************************** 

C  SUBRCUT 1NE  FCN2 
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C* ********  ****  *  *******  ****  ************  **************  ***** 
C 

SUBRCUT  INE  FCN2  (N  ,T,X,XPRIME1 
INTEGER  N 

REAL  X  (N1  •  XPR  IME(  M  »  T 
COMMON  /PARAM/  A ,  C  ,U  .  V  ,S  .  R  ,  B  ,  0 
XPRIME(11=-X( il*( L*A*X(21 1+R-B*X(2J 
XFRIME  (21  =-X(21*(  V*C*X  (11  1-D*X(11 
RETURN 
ENC 
C 

C*****  ****  *****4***  *?*  ************************  *  *  ***  *  ***** 
C  S  U ERG U 1  INE  FCN3 

C************* *************************************  ****** 

c 

SUBROUTINE  FCN3  (N  ,  T,  X,  XPR  IME1 
INTEGER  N 

REAL  >  (N1  »  XPR  IME(  N  J  »  T 
COMMON  /PARAM/  A.  C  ,U  ,  V  ,S  »  P  »B  ,  C 
XPRIME(1J=-X( 1 A  *i  L*A*X (21 I-B*X(2J 
XPRIME(2I=-X(2I*( V*C*X(1I 1-D*X(11 
RETURN 
ENC 
C 

£*************  ********  4- ***********************  *********** 

C  SLEROUINE  FCN4 

Q*  **** ********  *************************************  ****** 

c 

SUBROUTINE  FC N4 (N , T, X, XPR IME 1 
INTEGER  N 

REAL  MNJf  UPKIMEi  N)»T 
COMMON  /PARAM/  A,  C  ,U  ,  V  »S  ,  R ,  B  ,  C 
XPRIME(11--X(11*(U+A*X(21  I  *R—  B*X  (  2  I 
X  PRIME  (2 1  =  -X(  21*(  V*C*X(11  )+S-C*X( II 
RETURN 
ENC 
C 

£***** * *******  *  ***************************4***4********** 

C  SUBROUTINE  LP 

C* ****************************************** ************* 


729 

730 
405 


SUBROUTINE  LP  (  ALP  ,  EL  P  ,  CLP  .  SMM  IM 

INTEGER  I  ALP,  NLP,  Ml,  M2,IW  (631  »  I  EL P 

REAL  ALP ( 17,3  21  ,SLP, PSGL( 15! ,CSGL (151  f  JWI3E3J , 

SMMIN,YAL,ELP  (  151  ,  CLP  (151 

N  L  F=  1  5 

M  1=1  5 

M  2=0 

I  A LP=  1 7 

CALL  Z  >  A  LP  (AL  P  i  I A  L  P ,  B  LP » C  LP ,  N  LP  ,M  1  ,  M2  ,  SLP  ,  f  SOL 

fCSQLjRW,IW»IELPl 

VAL=(  l./SLPJ  -SMM1  N 

WRITE  (  1.72S1VAL 

FORMAT  (*0*  ,*  VALUE  OF  GAME-='  ,  F  14.3  1 
DC  405  JJ»1, 1  5 

P S C L  ( J J 1=P SOLI JJ1/SLP 
DSC  L ( J  J 1=0  SOL ( J Jl /SLP 
WRITE(  1  ,7301 JJ ,FSCL(JJ1»JJ 
FGFMAT  (  *0'  »  2X,  •  FRCe  OF.STR 
,F  12.3,'PROB.GF  STR.  #'»I3 
CONTINUE 
RETURN 
ENC 


,CSGL ( JJ  1 

*  ',13,'OF  Y=' 
•UF  X=*  ,F12. 31 


121 


APPENDIX  I 


TRANSFORMING  MIXED  STRATEGY  PROBLEM  INTO  LINEAR  PROGRAMMING 


The  payoff  function  has  been  defined  to  be 

A(X , Y)  =  Lv  -  Lx 

Y  selects  his  optimum  mixed  strategies  which  yeild 


rr  1 

V  ==  ) max 

qj 


“  m  m  m 

. 

-  J-l  J=l  J=) 


where 


payoff  to  x  when  x  adopts  ith  strategy  and  y  adopts 
j th  strategy 

probability  that  y  selects  j th  strategy 


Subject  to  : 


Bj- 


.=1,  q •  >  0  ,  j  =  1,2, 


Let  Vg  -  max 


. 

Lj=1  JM  j*t  J 


Then  the  original  problem  becomes 


minimize 


Subject  to  : 


Subject  to  : 


>  1 
j 

DVj  » 1 


S3  mj 

j 


>  1 


Qj  >  o  ,  j 


1  qi 

Since  Q0  =  -  and  Q.  =  -f- 


0 


J  v 


0 


~  1 , 2 ,  .  .  .  ,m 


””  qi  ‘  9iv°  =  % 


After  solving  the  LP  problem,  the  optimum  strategies  for  y 

*  * 

is  given  by  Vq  Some  constants,  K  =  |min(a^)|  could 

have  been  added  to  a^  to  ensure  Vq  >  0.  If  this  is  done, 

K  has  to  be  subtracted  from  the  optimum  value  obtained  by  the 
LP,  that  is 


* 

v 


K 


* 

where  v  =  the  value  of  the  game. 
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APPENDIX  J 


BACKGROUND  DATA  ON  KOREAN  WAR 


PERIOD  CONSIDERED  :  25  June  1950  to  7  July  1950 

(No  American  ground  involvement 
as  yet) 

TYPE  OF  ENGAGEMENT  :  Predominantly  land  combat 
GENERAL  STATE  OF  READINESS  : 

NORTH  KOREA  :  Well  prepared  by  1950;  arms 

build-up  and  training  of 
troops  since  1945;  many 
military  leaders  and  combat 
personnel  were  war  veterans 
fighting  in  China 

REPUBLIC  OF  KOREA  :  By  1950;  a  small  defense 

force  began  to  take  shape 
through  American  aid; 
training  only  started  around 
1948. 


SOURCE  OF  DATA  :  a)  Appleman,  R.E.,  United  States 

Army  in  the  Korean  War,  Department 
o£  the  Army,  1961 . 

b)  Montross,  Lynn,  U.S.  Marine  Opera- 
tions  in  Korea,  0  .~S~.  Marine  Corps.  , 
1954  . 

RELATIVE  STRENGTH  : 

NORTH  KOREA  REPUBLIC  OF  KOREA 


Total  strength  =  Q  =  135,000  men  Q  =  95,000  men 

x  y 


Tanks  :  150 

Artillery  pieces  :  1,600 


nil 

700 


Aircraft 

(i)  fighters  -  40 

(ii)  attack  bombers  -  70 

(iii)  reconnaissance  -  10 


no  combat  aircraft 
(22  trainer,  4 
auxiliary;  no  pilot) 


6 


CORRESPONDING  PARAMETERS  USED  IN  MODEL 


IIST  Of  EEFERENC ES 


Dolansky.  L.,  Present  State  of  the  Lancaster  Theory 
cf  Combat .  0 peratIons"Hesearch  Society  or  A  mer lea ,  T2. 


Rozencraf t ,  J. M.and  Moose,  P.H..  Lanchaster.  Theory  an 
Matrix  Games,  froc.  of  the  sixth  HTT7Cni“^orxsTloy  on 
sysleis7”15ecember  198  3. 


Siliak.  D.D.,  Non-linear  Systems,  John  Wiley  &  Sons, 
1569.  - - 


Ogata,  K.,  Modern  Control  Engineering,  Prentice-Hall, 
1970. 


Taylor,  James  G.,  Lanchester  Models  of  War  fare. 
Operations  Research  Society  ol'ffmerlca,  19B3. 


Lanchester,  F.R.,  Aircraft  in  Warfare  :  The  Down  of 
Fourth  Arm  No.  V.,  "lEe"  principle  ol  Concentration  , 
Engineer In g“S 0 ,  U22-4  23,  1914. 


Strang,  1. , 
Academe  Press, 


linear  Algebra  and 
T9757- pa ges-yi  9 -3T5 


Its  A pplicaticns 


Richter,  Stephen  L.,  Continuation  Methods  :  Theory  and 
Applications,  IEEE  "Trans.  on  system,  Kan  ANT5 
Cyrer nefics",  vol.  SMV-13,  No.  4,  July/August  198  3. 


Ohihiro,  H.,  Non  Linear  Oscillations  i n  Fhy sical 
System.  Mcgrav  Hill,  T951." 


Taha.  H.A.,  Cperati  ons  Research  :  An  Introduction. 
MacMillan,  1971. 


Frederick,  S.H.  and  Gerald,  J.L.,  Opera tion s  Research. 
Holden-Day,  1974. 


Liu,  C.L..  Introduction  to  Combinatorial  Mathematics, 
Mcgraw-Hill,~TTEB.  "  " 
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